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INTRODUCTION 


Magnetic  resonance  imaging  (MRI)  offers  a  powerful,  non-invasive  means  to  image  the  human  body.  Recent 
advances  in  imaging  the  brain  have  led  to  the  potential  of  monitoring  iron  build  up,  blood  products  in  dementia, 
stroke,  traumatic  brain  injury,  and  the  measurement  of  oxygen  saturation.  These  advances  relate  directly  to 
both  monitoring  and  measuring  the  amount  of  iron  in  brain  tissues  that  have  changes  in  magnetic  susceptibility 
of  each  individual  tissue.  Therefore,  the  proper  quantification  of  iron  will  play  a  key  role  in  future  radiology 
practice.  While  a  number  of  groups  using  several  methods  are  approaching  this  topic,  the  current  problems  are 
that  no  standard  protocol  has  been  established  nor  sensitivities  and  uncertainties  of  these  methods  have  been 
properly  evaluated.  These  methods  include  the  background  phase  removal  techniques.  Without  the  standard 
protocol  or  understanding  of  uncertainties  of  developed  methods,  the  quantified  susceptibility  values  can  be 
questionable.  As  a  first  step,  we  want  to  establish  a  standard  protocol  with  accurate  methods  for  the  magnetic 
susceptibility  quantifications  of  materials  in  phantoms.  We  have  accurately  quantified  susceptibilities  of  several 
concentrations  of  four  different  materials  in  water  with  MRI.  The  detailed  methods  and  results  are  presented 
below. 

BODY 

I.  Phantom  constructions 

Our  first  step  is  to  construct  gel  phantoms  inserted  with  straws,  in  order  to  establish  standards  in  the  magnetic 
susceptibility  quantification.  Currently,  we  know  that  we  can  apply  the  2D  CISSCO  (Complex  Image 
Summation  around  a  Spherical  or  a  Cylindrical  Object)  method  [Cheng  et  al.,  2009],  direct  phase 
measurements  (e.g.,  [Shen  et  al.,  2008]),  or  R2’  measurements  (e.g.,  [Haacke  et  al.,  2005])  for  susceptibility 
quantification  of  a  given  cylindrical  object.  In  order  to  avoid  partial  volume  and  Gibbs  ringing  artifacts,  the  direct 
phase  and  R2’  measurements  are  better  to  be  applied  to  cylindrical  objects  with  diameters  of  more  than  6 
pixels.  The  2D  CISSCO  method  breaks  through  that  limit. 

We  have  obtained  clear  straws  with  three  different  diameters  (roughly  3  mm,  5  mm,  and  6  mm)  and  we  also 
have  straws  with  4  different  colors  but  all  with  diameters  of  4  mm  from  the  same  manufacturer.  It  turned  out 
that  it  was  not  that  easy  to  obtain  clear  straws  with  different  diameters  from  the  same  manufacturer,  but  we 
tried  to  do  that  in  order  to  satisfy  reviewers’  request  for  this  project.  A  general  worry  is  that  some  straws  may 
contain  substances  to  influence  our  susceptibility  quantification. 

For  our  gel  phantoms,  we  choose  agar  as  the  material.  There  are  at  least  5  different  materials  to  make  gel. 
Agar  has  one  advantage  as  it  becomes  solid  during  the  room  temperature  but  gelatin  (another  material)  will 
melt.  However,  this  also  means  that  we  need  to  work  with  agar  at  a  relatively  high  temperature  (above  85  °C) 
and  straws  can  melt  at  that  temperature.  Our  experiments  show  that  we  can  insert  straws  with  liquid  samples 
into  agar  solutions  at  around  75  °C  without  melting  the  straws.  We  also  need  to  mix  agar  well  when  it  is  in  the 
solution  form.  We  have  purchased  a  4”  long  magnetic  stir  bar  and  it  serves  the  purposes  well. 

For  straws  that  are  narrower  than  5  mm,  it  is  also  a  problem  to  place  water  or  gel  solution  into  those  straws, 
due  to  the  surface  tension.  We  overcome  this  problem  by  covering  both  ends  of  each  straw  with  wax.  We  pick 
a  hole  through  wax  at  one  of  the  ends  and  inject  the  liquid  sample  from  the  bottom  end  of  each  straw.  Then  we 
seal  the  hole  with  cement. 

We  have  imaged  several  phantoms  in  order  to  examine  the  susceptibility  effect  due  to  the  walls  of  straws.  We 
fill  each  straw  with  the  distilled  water,  as  we  expect  minimal  susceptibility  difference  between  the  distilled  water 
and  agar  gel.  In  the  following  images  and  their  associated  tables,  we  show  results  from  three  phantoms 
acquired  from  the  shortest  and  the  longest  echo  times  of  an  11 -echo  SWI  (susceptibility  weighted  imaging) 
sequence  from  our  3T  MRI  machine.  The  imaging  resolution  was  1  mm  x  1  mm  x  1  mm.  The  in-plane  matrix 
size  was  256  x  256.  A  total  of  72  slices  were  acquired  but  only  64  were  reconstructed,  due  to  the  choice  of 
12.5%  oversampling.  Repetition  time  TR  was  37  ms  and  echo  times  ranged  from  5.68  ms  to  30.51  ms.  Echo 
spacing  was  about  2.5  ms.  Flip  angle  was  15°.  The  scan  time  was  11  min  and  22  sec.  A  32  x  32  high  pass 
filter  has  been  applied  to  phase  images  in  the  following  examples.  At  the  longest  echo  time,  we  can  see  very 
minor  susceptibility  effects  from  some  straws,  due  to  the  walls  of  the  straws.  Some  quick  estimations  of  the 
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magnetic  moments  of  the  walls  reveal  that  the  values  are  around  3-5  radian*pixel2.  If  we  carefully  examine  the 
phase  images,  in  general  we  see  that  smaller  straws  show  some  phase  patterns  in  images.  Thus,  we  choose 
to  use  straws  with  diameters  of  roughly  6.3  mm  for  our  subsequent  experiments,  as  those  straws  do  not  show 
any  obvious  susceptibility  effect  from  the  walls  of  straws. 


Phantom  1:  (Slice  no.  33) 
Straw  info: 


Straw  No. 

1 

2 

3 

Color 

Clear 

Green 

Orange 

Diameter  (mm) 

5.99  ±  0.06 

3.94  ±  0.05 

3.90  ±  0.05 

Echo  time:  5.68  ms 


(a)  magnitude 


(b) phase 
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Phantom  2:  (Slice  no.  33) 
Straw  info: 


Straw  No. 

1 

2 

3 

Color 

Clear 

Pink 

Green 

Diameter(mm) 

2.63  ±  0.03 

3.96  ±0.05 

3.94  ±  0.05 

Echo  time:  5.68  ms 


(e)  magnitude 


(f)  phase 


Echo  time:  30.51  ms 


(g)  magnitude  (h)  phase 
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Phantom  3:  (Slice  no.  32) 
Straw  info: 


Straw  No. 

1 

2 

3 

Color 

Pink 

Yellow 

Clear 

Diameter(mm) 

3.96  ±  0.05 

3.91  ±0.05 

2.63  ±  0.03 

Echo  time:  5.68  ms 


(i)  magnitude  (j)  phase 

Echo  time:  30.51  ms 


(k)  magnitude  (I)  phase 
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II.  Susceptibility  and  relaxation  rate  measurements  from  MRI  and  concentration  measurements  from 
ICP-AES 

1)  Sample  preparations 

We  prepared  several  concentrations  of  nanoparticle,  gadolinium,  calcium,  and  ferritin  solutions  for 
measurements.  Based  on  our  experiences  and  results  in  the  literature  [Shen  et  al.,  2008;  Weisskoff  and 
Kiihne,  1992;  Zheng  et  al.,  2013]  or  from  sigma.com,  the  highest  concentrations  of  nanoparticle,  gadolinium, 
calcium,  and  ferritin  solutions  prepared  by  us  were  aimed  at  25  mg/L  of  iron  in  nanoparticles,  927  mg/L  of  MR 
gadolinium  contrast  agent,  525  g/L  of  calcium,  and  1455  mg/L  of  iron  in  ferritin,  respectively.  Except  for 
nanoparticle  solutions,  these  values  would  correspond  to  +2  ppm  of  magnetic  susceptibility  values  relative  to 
the  susceptibility  of  water  (which  is  -9.05  ppm  in  SI  system).  For  nanoparticles,  the  25  mg/L  of  iron  would 
correspond  to  a  relative  susceptibility  value  of  +1  ppm.  In  order  to  achieve  these  expected  values,  we  first 
diluted  the  originally  very  high-concentrated  nanoparticle  solution,  the  MR  contrast  agent,  and  ferritin  solution 
with  distilled  water.  After  creating  those  large-quantity  batches,  we  further  repeatedly  diluted  each  solution  with 
the  distilled  water  by  half  of  its  previous  concentration  such  that  we  produced  at  least  four  different 
concentrations  for  each  of  the  four  materials.  For  MRI  experiments,  we  injected  sufficient  volume  of  each 
sample  into  a  clear  straw  with  a  diameter  of  roughly  6.3  mm  and  sealed  both  ends  of  the  straw.  For  each 
material,  we  placed  4  straws  with  different  concentrations  in  one  agar  phantom.  The  remaining  amount  of  each 
sample  was  stored  in  test  tubes  for  inductively  coupled  plasma  atomic  emission  spectrometer  (ICP-AES), 
magnetometer,  and  any  other  future  studies  (when  needed).  Our  plan  was  to  measure  each  sample  by  MRI, 
ICP-AES,  and  a  SQUID  (superconducting  quantum  interference  device)  based  magnetometer. 

2)  ICP-AES  measurements 

As  our  solutions  have  metallic  concentrations  in  or  above  the  mg/L  range,  it  is  more  appropriate  to  measure 
our  samples  with  ICP-AES.  ICP-AES  requires  2  ml  of  each  sample  and  uses  standard  metallic  concentrations 
as  references  to  establish  calibrations.  Even  at  a  very  high  concentration  such  as  850  mg/L,  the  quantified 
concentration  can  be  differed  by  1%  from  time  to  time.  When  the  concentration  is  between  10  mg/L  and  1 
mg/L,  the  uncertainty  is  at  least  10%.  Below  1  mg/L,  the  uncertainty  is  more  than  100%.  We  originally 
proposed  to  use  inductively  coupled  plasma  mass  spectrometer  (ICPMS)  to  measure  metallic  concentrations  in 
solutions.  However,  ICPMS  requires  the  concentration  of  each  metallic  content  to  be  less  than  0.2  mg/L.  This 
means  that  we  must  dilute  each  sample  for  ICPMS  measurements.  This  does  not  seem  a  problem  in  theory, 
but  we  prefer  (and  almost  insist)  to  measure  exactly  the  same  sample  by  different  techniques.  In  fact,  we  later 
found  out  from  ICP-AES  measurements  that  when  the  dilution  of  a  sample  was  by  a  factor  of  100,  we  easily 
observed  at  least  10%  deviation  from  the  concentration  of  the  undiluted  sample. 

For  nanoparticle  measurements  by  ICP-AES,  we  were  told  that  we  need  to  "digest"  the  samples  i.e.,  to  break 
iron  atoms  inside  each  nanoparticle.  Otherwise  the  measured  concentration  can  be  far  less  than  the  actual 
concentration  of  iron.  This  means,  we  need  to  add  concentrated  nitric  acid  to  each  of  our  nanoparticle  samples 
for  hours  or  days,  and  then  dilute  each  sample  such  that  the  concentration  of  nitric  acid  is  less  than  5%.  We 
have  measured  several  nanoparticles  concentrations  by  ICP-AES  with  and  without  “digestion.”  Surprisingly, 
ICP-AES  consistently  showed  that  each  digested  solution  had  a  lower  concentration  than  its  undigested 
counterpart.  This  is  completely  contradicting  to  theory  and  we  suspect  that  the  problem  may  be  due  to 
dilutions.  In  the  results  below,  we  show  the  measurements  of  ICP-AES  from  undigested  nanoparticle  solutions. 
Digestion  is  not  needed  for  the  gadolinium  MRI  contrast  agent  or  ferritin.  However,  as  the  concentrations  of 
ferritin  solutions  are  higher  than  the  ICP-AES  detection  limits  and  the  ICP-AES  machine  is  currently  under 
repair,  we  will  simply  rely  on  the  concentration  of  the  reference  ferritin  solution  purchased  from  sigma.com. 
Similarly,  the  concentrations  of  calcium  solutions  are  very  very  high.  We  thus  can  estimate  the  concentrations 
by  weighing  the  amount  of  CaCI2  (also  purchased  from  sigma.com)  and  the  amount  of  distilled  water.  The 
uncertainty  of  each  estimated  concentration  for  calcium  solutions  is  about  10%. 

3)  MRI  parameters  and  considerations  of  image  analysis 

We  have  imaged  our  phantoms  on  a  3  T  Siemens  Verio  machine  with  a  3D  11 -echo  SWI  sequence  and  a  2D 
14-echo  spin  echo  sequence  without  an  echo  train  length.  These  two  sequences  allowed  us  to  quantify 
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susceptibility  and  relaxation  rates  R2*  and  R2,  respectively.  As  a  12-channel  radiofrequency  (rf)  coil  was  used, 
it  would  be  better  to  analyze  images  from  each  rf  channel.  Our  best  choice  was  to  analyze  the  k-space  data  of 
four  channels  output  from  3  T  software.  We  also  imaged  each  phantom  with  these  sequences  at  two 
orientations,  when  straws  were  perpendicular  and  parallel  to  the  main  field.  The  in-plane  imaging  resolutions 
for  both  sequences  were  1  mm  x  1  mm.  The  in-plane  matrix  size  was  256  x  256.  For  SWI,  repetition  time  TR 
was  37  ms  and  echo  times  ranged  from  5.68  ms  to  29.48  ms,  with  an  echo  spacing  of  2.38  ms.  Flip  angle  was 
15°.  A  total  of  72  slices  with  slice  thickness  1  mm  were  acquired  but  only  64  were  reconstructed,  due  to  the 
choice  of  12.5%  oversampling.  The  scan  time  was  11  min  and  22  sec.  A  32  x  32  high  pass  filter  with  the 
Hanning  window  to  smooth  data  in  k-space  was  applied  to  phase  images  as  one  of  the  means  to  remove  the 
background  phase  before  susceptibility  quantifications.  Later  we  will  show  results  from  other  background 
phase  removal  method.  For  the  14-echo  spin  echo  sequence,  the  echo  time  ranged  from  12  ms  to  168  ms  with 
an  echo  spacing  of  12  ms.  The  repetition  time  TR  was  300  ms.  Only  the  middle  4  slices  with  slice  thickness  of 
4  mm  were  acquired.  The  scan  time  per  slice  was  about  76  sec  such  that  the  total  scan  time  for  this  spin  echo 
sequence  was  about  5  min  4  sec.  We  have  noted  that  the  actual  field  strength  of  our  Verio  machine  is  2.89  T. 

To  quantify  susceptibility  values,  as  a  standard  procedure,  we  typically  applied  our  2D  CISSCO  method  to 
straws  from  the  central  slice  of  an  SWI  dataset.  We  applied  the  CISSCO  method  to  other  slices  when  we 
wanted  to  check  whether  the  quantified  susceptibility  values  were  consistent  throughout  the  entire  phantom. 
The  current  CISSCO  method  provides  us  the  measurements  of  magnetic  moment  and  its  uncertainty  from 
each  straw.  Then  we  can  calculate  the  susceptibility  value  of  each  sample  by  dividing  the  magnetic  moment  by 
cross  section  of  each  straw,  as  we  have  measured  the  inner  diameter  and  outer  diameter  of  each  type  of  straw 
with  a  Vernier  calibrator.  In  addition,  for  our  susceptibility  quantification  of  an  intended  sample  inside  a  straw, 
we  tried  to  choose  the  phase  image  from  an  rf  channel  with  the  highest  intensity  around  the  straw  in  its 
corresponding  magnitude  image.  This  is  to  reduce  the  Gaussian  noise  around  the  straw  of  interest  in  the 
phase  image.  However,  the  final  magnitude  image  for  the  CISSCO  analysis  was  the  sum-of-squares 
magnitude  image  from  all  four  channels.  This  is  to  minimize  the  magnitude  intensity  variation  around  the  straw 
and  the  error  of  susceptibility  quantification. 


4)  Susceptibility  quantification  of  Fe304  nanoparticles 


Figures  1  and  2  show  two  phantoms  with  4  straws  and  3  straws,  respectively,  filled  by  7  different 
concentrations  of  Fe304  nanoparticles.  The  high  pass  filter  has  been  applied  to  phase  images.  We  quantified 
the  susceptibility  values  of  straws  directly  from  complex  images  with  TE  of  15.24  ms  and  we  were  able  to 
achieve  reasonably  accurate  measurements,  as  listed  in  Table  1.  The  phase  pattern  of  the  5th  straw  in  Fig.  2 
also  indicates  that  we  can  quantify  its  susceptibility.  We  do  not  show  that  result  here  due  to  the  large  error.  If 
we  can  properly  remove  the  background  field  from  phase  images  at  longer  echo  times,  we  can  reduce  the 
uncertainty  of  each  measured  susceptibility  value  from  our  CISSCO  method.  Similarly,  the  uncertainty  of  a 
quantified  susceptibility  value  is  smaller  if  the  susceptibility  value  is  larger. 

. KSS5 


Figure  1:  Coronal  magnitude  (left  panel)  and  filtered  phase  (right  panel)  images  of  nanoparticles  with 
concentrations  #1 ,  #2,  #3,  and  #4.  This  set  of  image  was  acquired  from  echo  time  15.24  ms  and  slice  No.  33. 
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Figure  2:  Coronal  magnitude  (left  panel)  and  filtered  phase  (right  panel)  images  of  nanoparticles  with 
concentration  #5,  #6,  and  #7.  This  set  of  image  was  acquired  from  echo  time  15.24  ms  and  slice  No.  33. 

The  quantified  results  from  ICP-AES  and  MRI  for  the  first  4  straws  are  listed  in  Table  1 .  Given  our  experimental 
findings  mentioned  above,  we  assume  no  susceptibility  effect  from  the  straw  walls  for  the  susceptibility  (Ax) 
quantification.  The  linear  fit  between  the  ICP-AES  and  MRI  susceptibility  measurements  is  shown  in  Fig.  3.  For 
the  straws  shown  in  Fig.  2,  we  have  tried  to  select  different  echo  times  for  susceptibility  quantifications  of 
different  straws,  in  order  to  obtain  sufficient  phase  distributions  around  each  straw.  Reliable  quantification  of 
low  susceptibility  values  from  those  straws  requires  a  background  removal  procedure  that  we  continue  to  work 
on. 


No. 

Cone  (mg/L) 

Ax  (ppm) 

S(Ax)/AX  (%) 

1 

25.70 

1.08 

3.8 

2 

12.59 

0.56 

4.3 

3 

6.80 

0.25 

8.3 

4 

3.00 

0.12 

17.3 

Table  1:  Quantified  results  from  ICP-AES  and  MRI  for  different  concentrations  of  nanoparticles.  The  first 
column  labels  each  straw.  The  second  column  lists  the  concentrations  of  iron  atoms  measured  by  ICP-AES. 
These  concentrations  are  different  from  the  concentrations  of  nanoparticles  (which  is  Fe304)  but  can  be 
converted.  The  third  column  Ax  lists  the  quantified  susceptibility  values  from  CISSCO.  The  fourth  column  lists 
the  percentage  error  of  each  susceptibility  measurement  derived  from  the  CISSCO  method. 

As  the  CISSCO  method  currently  only  uses  phase  distributions  outside  each  straw  for  the  susceptibility 
quantification,  in  order  to  check  consistency  of  our  results,  we  also  checked  whether  the  phase  value  inside 
each  straw  agree  with  that  calculated  from  the  quantified  susceptibility.  From  both  orientations  (parallel  and 
perpendicular  to  the  main  field),  we  have  discovered  that  phase  values  inside  many  straws  were  not  well 
predicted  by  the  quantified  susceptibility  values.  Closer  examines  revealed  that  phase  values  inside  different 
straws  in  Fig.  1  appeared  to  be  somewhat  proportional  to  the  quantified  susceptibility  values.  Similar  situation 
was  also  observed  from  the  parallel  orientation.  However,  phase  values  inside  straws  from  the  parallel  case 
were  much  less  than  phase  values  from  the  perpendicular  case  for  those  nanoparticle  solutions.  This  is 
completely  opposite  to  theoretical  predictions.  We  had  suspected  materials  inside  straws  other  than 
nanoparticles,  as  we  recalled  that  those  nanoparticles  were  mixed  with  sodium  hydroxide  (NaOH)  for  coating 
purposes  during  their  manufacture  processes.  However,  when  we  imaged  the  NaOH  solution  inside  a  straw, 
we  did  not  observe  any  substantial  phase  value  inside  the  straw,  when  the  straw  was  perpendicular  to  the 
main  field.  Thus,  we  have  carefully  examined  this  issue  by  simulating  the  effect  due  to  the  32  x  32  high  pass 
filter  applied  to  the  nanoaprticle  concentrations,  locations  of  straws,  and  phantom  sizes.  In  addition,  we  have 
examined  the  outcomes  due  to  different  background  phase  removal  methods.  We  will  present  those  results 
later  in  this  report. 
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Figure  3:  Quantified  susceptibility  values  from  MRI  versus  concentrations  of  nanoparticles  from  ICP-AES. 

5)  Susceptibility  quantification  of  gadolinium  Magnevist  MRI  contrast  agent 

Similar  to  nanoparticle  experiments,  we  applied  the  same  procedures  to  gadolinium.  Examples  are  shown  in 
Figs.  4  and  5  and  results  are  listed  in  Table  2  and  plotted  in  Fig.  6.  We  again  see  a  good  linear  fit  between  the 
susceptibility  and  concentration  measurements.  We  also  note  that  the  susceptibility  and  concentration  is 
convertible,  as  the  atomic  weight  of  gadolinium  is  157.25  g  (per  mole)  and  the  molar  susceptibility  for 
gadolinium  is  0.339  ppm/mM  provided  by  Weisskoff  and  Kiihne  [1992],  With  the  conversion  factor  calculated 
from  these  numbers,  our  susceptibility  measurements  agree  very  well  with  the  concentrations  measured  from 
ICP-AES  listed  in  Table  2. 


Figure  4:  Coronal  magnitude  (left  panel)  and  filtered  phase  (right  panel)  images  of  gadolinium  with 
concentrations  #1,  #2,  #3,  and  #4.  This  set  of  image  was  acquired  from  echo  time  15.24  ms  and  slice  No.  33. 
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Figure  5:  Coronal  magnitude  (left  panel)  and  filtered  phase  (right  panel)  images  of  gadolinium  with 
concentration  #5,  #6,  and  #7.  This  set  of  image  was  acquired  from  echo  time  15.24  ms  and  slice  No.  33. 


No. 

Cone  (mg/L) 

Ax  (ppm) 

8(Ax)/AX  (%) 

1 

855.03 

1.78 

3.6 

2 

467.59 

0.95 

3.7 

3 

213.62 

0.52 

4.6 

4 

109.03 

0.28 

8.0 

Table  2:  Quantified  results  from  ICP-AES  and  MRI  for  different  concentrations  of  gadolinium  MRI  contrast 
agent.  The  first  column  labels  each  straw.  The  second  column  lists  the  concentrations  of  gadolinium  atoms 
measured  by  ICP-AES.  The  third  column  Ax  lists  the  quantified  susceptibility  values  from  CISSCO.  The  fourth 
column  lists  the  percentage  error  of  each  susceptibility  measured  from  the  CISSCO  method. 


Figure  6:  Quantified  susceptibility  values  from  MRI  versus  concentrations  of  Gd-DTPA  from  ICP-AES. 
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6)  Measurements  of  relaxation  rates  for  nanoparticle  and  gadolinium  solutions 

The  relaxation  rates  R2*  and  R2  are  calculated  from  the  11 -echo  SWI  and  the  14-echo  spin  echo  images, 
respectively,  with  an  in-house  software  SPIN.  From  the  statistical  point  of  view,  it  is  better  to  use  at  least  6 
echoes  for  R2*  or  R2  estimations.  However,  if  R2  or  R2*  is  large,  then  it  is  also  important  to  check  whether  the 
signal-to-noise  ratio  (SNR)  is  at  least  3:1  in  images.  As  R2  or  R2*  quantifications  rely  on  exponential  fit  to  the 
MR  signals,  it  does  not  make  sense  to  include  images  with  very  low  SNRs  for  analyses.  With  the 
concentrations  or  susceptibility  values  we  have  considered,  we  may  use  less  than  6  echoes  to  quantify  R2  or 
R2*.  In  the  phantoms  shown  in  Figs.  1  and  3,  we  only  use  1  pixel  inside  each  straw  for  R2  or  R2*  calculations. 
This  is  to  avoid  partial  volume  or  Gibbs  ringing  effect.  In  addition,  with  the  given  SNR  in  images  and  the 
number  of  pixels  for  measurements,  we  will  also  be  able  to  estimate  the  uncertainty  for  each  measurement 
through  error  propagations. 

In  addition  to  measurements,  we  are  also  interested  in  theoretical  predictions  of  R2  and  R2*.  Gillis  et  al.  [1999] 
and  Koenig  and  Kellar  [1995]  have  several  theoretical  formulas  for  the  estimations  of  R2  values  based  on 
susceptibility  values  or  concentrations.  However,  we  have  found  that  none  of  the  formulas  can  properly  predict 
the  R2  values  measured  from  our  nanoparticle  or  gadolinium  phantoms.  In  fact,  all  the  theoretical  values  are  at 
least  two  orders  of  magnitude  away  from  the  measured  R2  values.  Thus,  we  are  not  including  those 
confounding  results  in  this  report. 

On  the  other  hand,  by  neglecting  the  R2  value  of  water  or  gel,  we  may  be  able  to  estimate  theoretical  R2*theory 
value  given  by  [Yablonskiy  and  Haacke,  1994;  Shen  et  al.,  2008] 

R2*theory  =  0.4  YAXBo  (1) 

where  y  is  the  gyromagnetic  ratio  (271*42.58  MHz/T),  Ax  is  the  quantified  susceptibility  of  the  composite  sample 
or  material  (in  which  the  volume  fraction  has  been  included),  and  B0  is  the  main  field  strength  (2.89  T).  We 
tested  Eq.  1  with  nanoparticles  mixed  in  gel  in  Shen  et  al.  [2008].  With  those  straws  in  Figs.  1  and  3,  we  have 
quantified  their  R2  and  R2*  values  and  list  them  in  Tables  3  and  4. 


Cone  (mg/L) 

R2  (Hz) 

R2*  (Hz) 

R2  theory  (Hz) 

Straw  1 

25.70 

28.1  ±2 

128.2  ±11 

334.02 

Straw  2 

12.59 

22.2  ±  1.5 

109.3  ±9 

173.19 

Straw  3 

6.80 

14.7  ±2 

78.7  ±  8 

77.32 

Straw  4 

3 

11.7  ±2.5 

44.8  ±  7 

37.11 

Table  3:  Quantified  results  from  ICP-AES  and  MRI  relaxation  rates  with  uncertainties  for  different 
concentrations  of  nanoparticles  diluted  in  distilled  water.  The  first  column  labels  each  straw.  The  second 
column  lists  the  concentrations  of  iron  atoms  measured  by  ICP-AES.  The  third  column  lists  measured  R2 
values  from  the  spin  echo  sequence.  The  fourth  column  lists  R2*  values  measured  from  SWI.  These 
measurements  are  from  the  first  phantom  shown  in  Fig.  1.  Only  the  first  6  echoes  were  chosen  for  R2  and  R2* 
quantifications.  The  fifth  column  lists  the  theoretical  R2*  values  calculated  from  Eq.  1. 


Cone  (mg/L) 

R2  (Hz) 

R2*  (Hz) 

Straw  1 

855.03 

30.3  ±0.2 

31  ±2 

23.6 

Straw  2 

467.59 

15.4  ±0.1 

16  ±  2 

12.6 

Straw  3 

213.62 

5.7  ±0.7 

5.8  ±2 

6.9 

Straw  4 

109.03 

3.0  ±0.1 

3.2  ±2 

3.7 

Table  4:  Quantified  results  from  ICP-AES  and  MRI  relaxation  rates  with  uncertainties  for  different 
concentrations  of  gadolinium  diluted  in  distilled  water.  The  second  column  lists  the  concentrations  of  iron 
atoms  measured  by  ICP-AES.  The  fifth  column  lists  the  theoretical  R2  values  estimated  from  the  molar 
susceptibility  of  gadolinium,  quantified  susceptibility  values  (from  Table  2),  and  relaxivity  4.5  /mM/sec  from 
Haacke  et  al.,  Magnetic  Resonance  Imaging,  Ch.  22  [1999].  Other  columns  have  been  described  in  Table  3. 
Again,  only  the  first  6  echoes  were  chosen  to  quantify  R2  and  R2*  along  with  uncertainties  for  all  4  straws. 
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We  also  plot  the  results  as  a  function  of  susceptibility  values  (which  are  linearly  proportional  to  concentrations) 
and  show  them  in  Figs.  7  and  8.  Clearly  shown  in  both  Table  3  and  Fig.  7,  R2*theory  only  agrees  with  quantified 
R2*  of  low  concentration  nanoparticle  solutions  in  straws  #3  and  #4.  Between  the  listed  R2*theory  values  in 
Table  3  and  measured  R2*  values  in  Table  4,  it  is  also  clear  that  Eq.  1  cannot  predict  R2*  values  from 
gadolinium  solutions.  The  theoretical  R2  values  in  Table  4  are  based  on  the  empirically  measured  relaxivity  4.5 
/mM/sec  for  gadolinium  from  Haacke  et  al.,  Magnetic  Resonance  Imaging,  Ch.  22  [1999],  Given  the  estimated 
uncertainties  of  R2  values,  those  theoretical  R2  values  do  not  quite  agree  with  measured  values. 


Figure  7:  Quantified  R2  and  R2*  values  versus  quantified  susceptibilities  of  nanoparticles.  The  relation  between 
the  susceptibility  and  R2*  measurements  is  not  linear. 


Figure  8:  Quantified  R2  and  R2*  values  versus  quantified  susceptibilities  of  Gd-DTPA. 
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Nonetheless,  Fig.  8  shows  linear  relations  between  R2/R2*  and  susceptibility  values  (and  thus  concentrations) 
of  gadolinium  solutions.  All  these  results  prompted  us  to  further  investigate  the  R2  and  R2*  quantifications.  We 
have  prepared  a  few  more  phantoms  with  straws  and  tubes  briefly  described  below.  We  have  also  mixed 
nanoparticles  with  gelatin  solutions  and  have  made  solidified  samples. 

In  the  third  phantom,  we  prepared  a  straw  with  nanoparticle  solution  having  27.50  mg/L  of  iron  in  water  and 
two  plastic  tubes  of  diameter  14  mm  with  nanoparticles  mixed  in  the  gelatin  solutions.  For  tube  #1,  30  mg/mL 
of  gelatin  solution  was  prepared  and  mixed  with  27.50  mg/L  nanoparticle  solution  in  1:1  ratio  by  volume.  For 
tube  #2,  50  mg/mL  of  gelatin  solution  is  prepared  and  mixed  with  nanoparticles  solution  in  1:1  ratio  by  volume. 
The  concentrations  of  iron  nanoparticles  in  both  tubes  were  12.59  mg/L.  A  total  of  16  pixels  inside  each  tube 
were  used  for  R2  and  R2*  quantifications.  The  results  are  listed  in  Table  5.  We  now  see  the  agreement 
between  R2*theory  and  measured  R2*  for  solidified  samples.  However,  the  R2  or  R2*  values  from  those 
solidified  samples  are  statistically  significantly  higher  than  those  in  Table  3  with  the  same  concentration. 


Cone  (mg/L) 

R2  (Hz) 

R2*  (Hz) 

LU— Llil 

Straw  1 

25.70 

13.6  ±0.8 

1 12  ±  6 

334.02 

Tube  1 

12.59 

34  ±3 

159  ±9 

173.19 

Tube  2 

12.59 

34  ±3 

172  ±8 

173.19 

Table  5:  Relaxation  rates  of  nanoparticles  diluted  with  water  in  a  straw  and  mixed  with  gelatin  in  tubes. 
Nanoparticles  mixed  with  gel  were  solidified  in  tubes.  These  results  were  from  the  third  phantom.  For  R2  and 
R2*  quantifications,  the  first  6  echoes  were  used  for  the  straw  but  only  3  echoes  were  used  for  both  tubes,  due 
to  the  insufficient  SNR  after  the  third  echo.  The  meaning  of  each  column  was  given  in  Table  3. 

In  the  fourth  phantom,  the  three  straws  had  the  same  mixtures  as  in  the  straw  and  the  two  tubes  in  the  third 
phantom,  respectively.  The  results  are  listed  in  Table  6.  Again,  R2*the0ry  and  measured  R2*  values  agree  for 
solidified  samples  in  straws  #2  and  #3.  This  indicates  that  measurements  of  relaxation  rates  do  not  seem  to 
depend  on  geometry.  However,  there  may  be  a  minor  concern  about  different  R2*  values  between  Table  5  and 
Table  6  from  the  highest  concentration,  but  we  cannot  make  a  firm  conclusion  as  two  standard  deviations  of 
estimated  uncertainties  still  overlap  with  each  other.  Nonetheless,  the  R2  values  seem  to  be  affected  by  the 
gelatin  concentrations  but  this  is  inconsistent  between  straws  #2  and  #3  shown  in  Table  6.  In  addition,  we  do 
not  place  the  highest  nanoparticle  concentration  in  gelatin  form,  as  the  theoretical  R2*  value  (334  Hz)  indicates 
that  we  will  need  echo  times  between  1  ms  and  6  ms  with  echo  spacing  of  roughly  1  ms.  There  is  no  such  a 
practical  sequence  to  use. 


Cone  (mg/L) 

R2  (Hz) 

R2*  (Hz) 

R2  theory  (Hz) 

Straw  1 

25.70 

18  ±  1.2 

132  ±6.5 

334.02 

Straw  2 

12.59 

39  ±5.8 

159  ±  13 

173.19 

Straw  3 

12.59 

58.8  ±6.6 

177  ±  14 

173.19 

Table  6:  Similar  to  Table  5,  relaxation  rates  of  nanoparticles  in  water  or  in  solidified  gelatin  are  shown.  For  all  3 
straws,  the  first  6  echoes  are  chosen  to  quantify  R2*  and  their  corresponding  uncertainties.  However,  for  R2 
quantification,  only  the  first  3  echoes  are  used  for  straws  #2  and  #3  but  the  first  6  echoes  are  used  for  straw  1 . 

We  constructed  the  fifth  phantom  with  similar  containers  used  in  the  third  phantom.  We  placed  855.03  mg/L  of 
gadolinium  diluted  in  water  in  a  straw  and  mixed  gadolinium  with  the  1:1  ratio  by  volume  with  30  mg/mL  and  50 
mg/mL  of  gelatin  solutions,  respectively,  in  tubes  #1  and  #2  with  diameters  of  11  mm.  The  concentrations  of 
gadolinium  in  both  tubes  were  469.59  mg/L.  Table  7  lists  the  results.  The  results  in  the  straw  are  consistent 
with  the  results  shown  in  Table  4.  However,  the  results  from  tubes  with  solidified  gel  indicate  higher  R2  or  R2* 
values  as  their  counterparts  in  Table  4.  As  in  theory,  relaxation  rates  are  inversely  proportional  to  the  diffusion 
constant,  it  is  expected  that  the  relaxation  rates  measured  from  solidified  samples  (whose  diffusion  constant  is 
lower)  are  higher  than  those  measured  from  liquid  samples  with  the  same  concentrations  (whose  diffusion 
constants  would  be  higher). 
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Cone  (mg/L) 

R2  (Hz) 

R2*  (Hz) 

Straw  1 

855.03 

29.5  ±  0.2 

31.3  ±  1.6 

Tube  1 

467.59 

19.8  ±0.04 

21.4  ±0.5 

Tube  2 

467.59 

18.5  ±0.03 

19.9  ±0.4 

Table  7:  Relaxation  rates  of  gadolinium  diluted  with  water  in  a  straw  and  mixed  with  gelatin  in  tubes. 
Gadolinium  mixed  with  gel  were  solidified  in  tubes.  These  results  were  from  the  fifth  phantom.  For  R2  and  R2* 
quantifications,  the  first  6  echoes  were  used  for  analyses. 

We  placed  the  same  solutions  from  the  straw  and  tube  #1  used  in  the  fifth  phantom  into  the  straws  #1  and  #2 
in  the  sixth  phantom,  respectively.  We  also  placed  the  same  solution  in  straw  #2  in  a  tube  for  the  sixth 
phantom.  The  results  are  shown  in  Table  8  and  they  are  consistent  with  results  shown  in  Table  7.  We  further 
rotated  and  imaged  the  sixth  phantom  and  quantified  the  relaxation  rates  inside  the  straws  and  tube  again.  The 
results  listed  in  Table  9  are  consistent  with  those  shown  in  Table  8. 


Cone  (mg/L) 

R2  (Hz) 

R2*  (Hz) 

Straw  1 

855.03 

29  ±  0.2 

31.4  ±  1.4 

Straw  2 

467.59 

20  ±0.1 

22  ±  1.3 

Tube  1 

467.59 

20  ±  0.2 

21  ±  1.4 

Table  8:  Relaxation  rates  of  gadolinium  solutions  from  the  sixth  phantom.  For  quantification  of  relaxation  rates, 
the  first  6  echoes  have  been  used. 


Cone  (mg/L) 

R2  (Hz) 

R2*  (Hz) 

Straw  1 

855.03 

29.3  ±0.3 

30  ±  1.5 

Straw  2 

467.59 

20  ±  0.2 

21.3  ±  1.8 

Tube  1 

467.59 

19.2  ±0.1 

21  ±  1.7 

Table  9:  Quantification  of  relaxation  rates  from  the  sixth  phantom  with  a  90°  rotation  such  that  before  and  after 
the  rotation,  straws  and  the  tube  are  always  perpendicular  to  the  main  field. 

From  the  above  results  above,  it  is  apparent  that  relaxation  rates  of  either  nanoparticle  or  gelatin  solutions  are 
higher  in  solidified  gel  than  in  water.  However,  the  ratios  of  increased  relaxation  rates  are  not  the  same  for 
nanoparticle  and  gadolinium  samples.  This  means  that  a  simple  diffusion  constant  is  not  the  only  factor  to 
affect  the  changes  of  relaxation  rates  from  liquid  samples  to  solidified  samples.  In  addition,  for  nanoparticles, 
R2  and  R2*  are  very  different  but  R2  and  R2*  are  about  the  same  for  gadolinium  samples.  Furthermore,  the 
R2*  theory  can  only  predict  the  values  of  nanoparticles  in  solidified  gels.  R2  theory  cannot  predict  any 
measured  R2  values.  However,  when  the  gadolinium  samples  have  the  same  concentration  as  that  of  the 
second  straw  listed  in  Table  2,  regardless  whether  the  samples  were  mixed  with  water  or  gelatin,  the 
susceptibility  measurements  from  CISSCO  agreed  very  well  with  the  value  (within  the  uncertainty)  shown  in 
Table  2.  Similarly,  susceptibility  values  of  nanoparticle  samples  are  basically  the  same  for  the  same 
concentration,  independent  of  the  material  used  for  sample  preparations.  Some  results  are  listed  in  Table  10. 
All  these  findings  are  important  knowledge  to  the  MRI  community. 


Cone  (mg/L) 

AX  (ppm) 

Straw  1 

25.70 

1.28 

Straw  2 

12.59 

0.66 

Straw  3 

12.59 

0.68 

Table  10:  Quantified  susceptibility  values  from  the  fourth  phantom,  Table  6,  at  TE  =  8.07  ms,  with  the  use  of 
the  high  pass  filter  to  remove  the  background  phase.  As  the  susceptibility  measured  from  the  first  straw 
(containing  nanoparticles  in  distilled  water)  is  higher  than  the  susceptibility  of  the  first  straw  shown  in  Table  1, 
this  indicates  that  we  should  scale  the  quantified  susceptibility  values  of  straws  #2  and  #3  in  this  table  before 
we  compare  these  results  to  the  measured  susceptibility  of  straw  #2  in  Table  1.  This  is  because  lower 
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concentrations  were  always  diluted  from  the  highest  concentration  in  a  phantom.  After  this  scaling,  the 
quantified  susceptibility  values  of  straws  #2  and  #3  in  this  table  agree  very  well  with  the  susceptibility  of  the 
second  straw  in  Table  1. 

7)  Four  methods  for  removing  the  background  phase 

As  mentioned  above,  the  susceptibility  values  quantified  from  our  CISSCO  method  did  not  fully  agree  with 
results  derived  from  phase  values  inside  each  straw.  The  question  was  whether  the  disagreements  were  due 
to  post  processing  techniques  or  physical  phenomena  such  as  susceptibility  tensor.  We  first  examined  the 
images  of  the  phantom  (e.g.,  Fig.  4)  containing  4  straws  filled  with  different  concentrations  of  Gd-DTPA  in 
distilled  water.  The  nanoparticle  samples  contained  NaOH  which  may  further  complicated  the  susceptibility 
quantification  and  thus  would  be  examined  later.  We  began  our  investigations  with  the  use  of  the  32  x  32  high- 
pass  filter  and  then  other  background  removal  techniques  including  SHARP  [Schweser  et  al.,  2011],  quadratic 
fitting  (i.e.,  second  order  polynomial  fit  to  the  background  phase),  and  our  proposed  algorithm.  These  were  a 
total  of  4  methods.  It  is  necessary  to  remove  the  unwanted  background  phase  before  any  susceptibility 
quantification  method  can  be  applied  to  images.  In  order  to  use  SHARP  and  quadratic  fitting  method,  we  had 
to  unwrap  phase  in  phase  images.  After  applying  any  of  these  background  removal  techniques,  we  directly 
measured  the  averaged  phase  values  inside  each  straw  (which  contained  a  solution  at  a  given  concentration) 
from  different  echo  times.  Different  background  removal  methods  also  allow  us  re-quantify  the  susceptibility 
value  using  CISSCO  again.  Then  we  compared  these  results  for  consistency  checks. 

a)  High-pass  filter  and  simulations  of  its  effect 

The  high-pass  filter  has  been  widely  used  due  to  its  fast  post-processing  capability.  To  examine  the  high-pass 
filter  effect  on  susceptibility  quantification  and  phase  measurements,  we  have  conducted  various  simulations. 
In  those  simulations,  we  set  up  an  initial  40962  matrix  and  properly  generate  a  straw  with  a  radius  of  3  pixels 
centered  at  a  final  matrix  of  256  x  256.  We  either  simulate  the  actual  diameter  (125  pixels)  of  the  phantom  or 
fill  the  entire  matrix  space  by  the  gel  phantom.  The  quantified  magnetic  moments  from  former  and  latter  cases 
are  labeled  by  p”  and  p’,  respectively,  In  Fig.  9.  We  assign  a  series  of  magnetic  moments  p  to  straws  and 
generate  phase  images  with  straws  perpendicular  and  parallel  to  the  main  field.  In  both  perpendicular  and 
parallel  orientations,  we  apply  a  32  x  32  high-pass  filter  to  those  simulated  images.  For  the  perpendicular 
orientation,  we  quantify  each  magnetic  moment  value  p  by  CISSCO.  The  results  are  shown  in  Fig.  9  and  are 
compared  to  the  originally  assigned  p  values. 


*-p 


10  20  30  40  50  60  70 

p  (radians*pixelsA2) 


Figure  9:  Comparisons  between  magnetic  moment  values  before  and  after  the  high-pass  filter,  p’  denotes  the 
results  with  the  entire  matrix  space  filled  by  the  gel  phantom  and  p”  denotes  the  results  with  a  sizable  gel 
phantom,  p  indicates  the  input  value  we  assigned  in  each  simulation. 
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Figure  9  indicates  that  the  high-pass  filter  affects  nonlinearly  to  the  magnetic  moment  quantification.  In 
particular,  the  quantified  magnetic  moments  in  general  are  smaller  than  the  actual  values,  except  for  the  case 
when  p  =10  radians-pixel2.  The  error  induced  by  the  high-pass  filter  effect  is  within  7%  in  the  range  of  p  shown 
in  Fig.  9. 

To  further  study  the  background  phase  effect  in  the  gel  phantom,  we  have  also  added  a  constant  background 
phase  of  2  radians  inside  the  gel  phantom.  The  simulated  phase  images  for  p  =  10,  40,  and  70  radians-pixel2 
are  shown  in  Fig.  10.  We  apply  the  high-pass  filter  and  quantify  the  magnetic  moment  with  CISSCO.  These 
results  are  almost  the  same  as  those  obtained  from  zero  background  phase,  labeled  by  p”in  Fig.  9. 


Figure  10:  Simulated  phase  images  with  a  constant  background  phase  of  2  radians  inside  the  gel  phantom  for 
p  =  10,  40,  and  70  radians-pixel2.  A  straw  is  simulated  at  the  center  of  the  phantom. 

In  addition  to  the  quantification  of  magnetic  moment,  we  have  also  examined  the  high-pass  filter  effect  on 
phase  values  inside  straws.  From  simulated  phase  images  in  both  the  perpendicular  and  parallel  orientations, 
we  compare  the  input  phase  values  with  those  after  the  application  of  the  high-pass  filter,  either  with  a  sizable 
phantom  size  or  with  the  entire  field-of  view  (FOV)  filled  by  the  phantom.  Results  are  shown  in  Fig.  11 
(perpendicular  orientation)  and  Fig.  12  (parallel  orientation). 


Figure  11:  Comparisons  between  phase  values  inside  straws  before  and  after  the  high-pass  filter  in  the 
perpendicular  orientations.  4>’  denotes  the  results  with  the  entire  matrix  space  filled  by  the  gel  phantom  and  <J>” 
denotes  the  results  with  a  sizable  gel  phantom.  <|>  indicates  the  input  value  calculated  from  the  given  p  value  in 
each  simulation. 
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Figures  11  and  12  show  that  the  high-pass  filter  reduces  the  phase  value  inside  each  straw.  This  is  expected. 
However,  they  show  that  the  size  of  the  gel  phantom  has  to  be  included  in  simulations  when  the  high-pass  filter 
is  applied.  Otherwise  we  will  observe  a  large  difference  (or  error)  in  phase  measurements  for  the  perpendicular 
orientation  (see  Fig.  11).  This  is  a  big  surprise  to  us.  On  the  other  hand,  when  we  use  CISSCO  to  quantify  the 
magnetic  moment  with  the  application  of  the  high-pass  filter,  as  shown  in  Fig.  9,  the  size  of  the  phantom  is  not 
necessary  to  be  included  in  simulations.  Although  the  high-pass  filter  has  also  affected  the  phase  distributions 
outside  each  straw,  these  results  indicate  that  our  CISSCO  method  is  quite  robust  for  the  quantification  of  the 
magnetic  moment,  which  utilizes  the  phase  distributions  outside  each  straw. 


Figure  12:  Comparisons  between  phase  values  inside  straws  before  and  after  the  high-pass  filter  in  the 
parallel  orientation.  4>’  denotes  the  results  with  the  entire  matrix  space  filled  by  the  gel  phantom  and  f’  denotes 
the  results  with  a  sizable  gel  phantom.  <|>  indicates  the  input  value  calculated  from  the  given  p  value  in  each 
simulation. 

As  simulated  phase  values  f  and  tj>”  do  not  deviate  much  from  their  corresponding  actual  phase  values  in  Fig. 
12  (but  it  is  not  the  case  in  Fig.  1 1 ),  these  suggest  that  the  high-pass  filter  affects  more  on  areas  around  straws 
where  induced  phase  distributions  are  stronger  for  the  perpendicular  orientation.  Effectively,  those  areas  may 
be  treated  as  part  of  a  straw  by  the  high-pass  filter.  As  we  know  that  the  high-pass  filter  affects  more  on  larger 
objects,  non-uniform  deviations  between  4>’  and  f’  shown  in  Figs.  1 1  and  12  may  be  explained  by  this  reason. 

When  we  compare  high-pass  filtered  phase  values  inside  straws  from  actual  phantom  images  with  4>”  in  Figs. 
11  and  12,  only  phase  values  from  the  first  two  echoes  for  the  perpendicular  orientation  agree  well  with 
predicted  <J>”.  Closer  investigations  reveal  that  the  leftover  background  phase  in  images  after  the  high-pass  filter 
is  not  close  to  zero  for  images  from  longer  echoes  or  the  parallel  orientation.  This  prompts  us  to  consider  other 
background  phase  removal  methods. 

Our  high-pass  filter  results  based  on  CISSCO  were  quantified  from  TE  =  8.07  ms,  10.46  ms,  12.85  ms,  and 
17.63  ms  for  the  1st,  2nd,  3rd,  and  4th  straw,  respectively.  These  echo  times  were  chosen  in  order  to  optimize  the 
application  of  the  high-pass  filter  for  relatively  accurate  quantification.  The  results  are  given  in  Table  11  below 
and  they  can  be  considered  as  reference  values.  However,  for  phase  measurements  inside  straws  listed  in 
Table  12,  we  do  not  treat  those  values  as  references.  Figure  1 1  demonstrates  the  obvious  reason. 

b)  SHARP 

SHARP  is  a  state-of-art  tool  for  removing  the  background  phase.  As  this  is  a  method  developed  by  Schweser 
et  al.,  [2011]  we  have  re-built  the  algorithm.  For  SHARP,  we  could  take  the  images  from  the  longest  echo  time 
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at  TE  =  29.58  ms  and  unwrap  the  phase  in  images  with  the  prelude  algorithm  in  FSL.  With  the  unwrapped 
phase  images,  we  utilized  SHARP  in  both  perpendicular  and  parallel  orientation. 

In  the  perpendicular  orientation,  we  quantified  the  effective  magnetic  moment  and  susceptibility  for  all  4  straws 
and  compared  them  with  the  ICP-AES  measurements,  results  from  the  high-pass  filter,  and  results  from  our 
own  algorithm  (described  below  with  choices  of  echo  times).  Susceptibility  or  phase  values  of  ICP-AES  were 
calculated  from  their  measured  concentrations.  The  results  are  given  in  Table  11.  Both  SHARP  and  our  own 
algorithm  were  performed  at  the  echo  time  of  29.58  ms.  The  echo  times  used  for  the  high-pass  (HP)  filter  and 
for  different  straws  were  given  above. 

Table  11  shows  that  different  methods  led  to  similar  results.  For  comparisons,  we  performed  SHARP  with  two 
different  radii,  8  and  18  pixels,  labeled  by  SHARP-8  and  SHARP-18,  respectively.  Between  these  two  choices 
in  SHARP,  the  results  of  SHARP-18  gave  closer  values  to  the  susceptibility  values  derived  from  ICP-AES 
measurements.  The  uncertainty  of  each  measurement  from  CISSCO  is  no  more  than  5%  and  thus  is  not  listed. 
The  uncertainty  of  each  ICP-AES  measurement  is  also  small  and  has  been  mentioned  above.  The  accuracy  of 
our  proposed  algorithm  depends  on  the  accuracy  of  the  results  from  the  high-pass  filter  used  at  shorter  echo 
times.  Different  methods  agreed  better  at  higher  concentrations.  We  will  continue  to  investigate  those  slight 
differences  between  SHARP  and  other  methods. 


Straw 

#1 

#2 

#3 

#4 

ICP-AES 

1.84 

1.01 

0.46 

0.24 

HP  filter 

1.78 

0.95 

0.52 

0.28 

SHARP-8 

1.75 

0.90 

0.47 

0.22 

SHARP-18 

1.82 

0.92 

0.47 

0.25 

Proposed  algorithm 

1.74 

0.94 

0.55 

0.30 

Table  11:  Susceptibility  quantifications  (in  ppm)  of  different  methods  for  background  phase  removal  from  the 
perpendicular  orientation.  Each  straw  is  labeled  in  the  first  row. 

For  phase  values  inside  straws  at  the  perpendicular  orientation,  the  application  of  SHARP  leads  to  consistent 
results  between  measured  phase  values  of  gadolinium  solutions  and  quantified  susceptibility  values  using 
CISSCO.  However,  some  differences  remain  for  nanoparticle  solutions. 

To  further  examine  the  consistency  among  our  quantified  results,  we  also  applied  SHARP  to  parallel 
orientation  and  measured  the  phase  values  inside  each  straw  from  four  different  echo  times.  The  results  from 
TE  =  8.07  ms  are  listed  in  Table  12.  Phase  values  from  ICP-AES  were  calculated  from  measured 
concentrations.  As  the  SNR  was  high  in  each  image,  the  uncertainty  of  each  average  phase  value  listed  in 
Table  12  is  small  and  thus  neglected. 


Straw 

#1 

#2 

#3 

#4 

ICP-AES 

-2.45 

2.10 

0.96 

0.50 

HP  filter 

-1.62 

1.57 

0.78 

0.43 

SHARP-8 

-2.15 

2.04 

1.10 

0.58 

Quadratic  fitting 

-2.22 

2.05 

1.06 

0.58 

Table  12:  Phase  measurements  inside  straws  in  units  of  radian  from  different  methods  at  the  parallel 
orientation  at  TE  =  8.07  ms.  Again,  each  straw  is  labeled  in  the  first  row. 

Table  12  shows  that  the  high-pass  filter  can  significantly  alter  the  phase  values  when  the  background  phase 
distributions  are  strong.  SHARP  and  the  quadratic  fitting  method  give  more  reasonable  results  than  those  from 
the  high-pass  filter.  In  this  study,  SHARP-8  and  SHARP-18  give  almost  the  same  outcomes  so  we  only  show 
results  from  one  of  them.  Comparing  results  between  Table  11  and  Table  12  and  using  ICP-AES  results  as 
references,  it  seems  that  results  from  our  proposed  algorithm  are  the  most  consistent  with  results  from  ICP- 
AES.  For  example,  the  susceptibility  values  from  our  proposed  algorithm  are  less  than  values  from  ICP-AES 
for  straws  #1  and  #2  in  Table  1 1  but  more  than  values  from  ICP-AES  for  straws  #3  and  #4.  These  trends  agree 
with  results  in  Table  12  from  SHARP-8  and  quadratic  fitting  method,  if  we  can  believe  that  the  quadratic  fitting 
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method  has  also  properly  removed  the  background  phase.  However,  following  the  same  logic,  results  from 
SHARP  are  not  self-consistent  between  Table  1 1  and  Table  12. 

c)  Quadratic  fitting 

From  original  phase  images  of  the  parallel  orientation,  phase  aliasing  (wrapping)  can  occur  along  the  radial 
direction.  This  is  due  to  the  background  phase  from  the  gel  phantom  itself.  As  the  straws  parallel  to  the  field  do 
not  produce  magnetic  fields  outside  the  straws,  we  fit  the  phase  in  the  agar  gel  part  without  straws  to  a 
quadratic  equation.  Again,  we  first  unwrap  phase  aliasing  in  original  phase  images,  and  then  apply  binary 
masks  to  exclude  the  straws  and  the  noise  region  outside  the  agar  gel  phantom.  Then,  from  the  masked  phase 
images  we  are  able  to  perform  the  quadratic  fitting  method  and  remove  the  background  phase  with  a  complex 
division  procedure.  The  resulting  image  is  shown  in  Fig.  13c.  We  then  measure  the  phase  inside  each  straw  at 
different  echo  times.  Results  scaled  by  echo  times  are  consistent  within  10%  differences.  Results  from  TE  = 
8.07  ms  are  listed  in  Table  12.  For  the  perpendicular  orientation,  as  straws  induce  phase  distributions  outside, 
it  would  be  less  accurate  to  apply  the  quadratic  fitting  method. 


Figure  13:  (a)  Unwrapped  phase  from  the  original  image,  (b)  Background  phase  from  the  quadratic  fitting 
method,  (c)  Resulting  phase  image  after  removing  the  background  phase  (b)  from  (a). 

d)  Our  proposed  algorithm 

A  simple  algorithm  described  in  the  proposal  was  also  applied  in  the  perpendicular  orientation,  in  order  to 
remove  the  background  phase  for  susceptibility  quantification.  The  procedures  are  described  below: 

1.  We  applied  the  32  x  32  high-pass  filter  on  phase  images  from  a  shorter  echo  time  of  the  11 -echo  SWI,  for 
example,  TE  =  15.24  ms,  and  then  quantified  magnetic  moment  for  each  straw; 

2.  From  the  quantified  magnetic  moment,  we  simulated  phase  distributions  inside  and  around  each  straw; 

3.  We  applied  complex  division  to  remove  the  simulated  phase  pattern  around  each  straw  (step  2)  from  the 
filtered  phase  in  step  1 .  This  was  to  obtain  phase  images  containing  only  “background  phase”; 

4.  We  multiplied  the  “background”  phase  pattern  by  an  integer,  e.g.  2,  to  obtain  a  “background”  phase  at  TE  = 
30.48  ms  (2*15.24  =  30.48); 

5.  We  complex  divided  the  "background"  phase  (step  4)  from  the  phase  image  at  a  longer  echo  time,  e.g.  TE  = 
29.58  ms,  and  then  applied  an  8  x  8  high-pass  filter. 

6.  We  obtained  a  phase  image  at  TE  =  29.58  ms,  with  most  of  the  background  phase  removed  but  phase 
induced  by  the  straws  remained.  We  quantified  the  magnetic  moment  for  each  straw,  calculated  the 
susceptibility,  and  compared  it  with  results  from  other  methods.  The  results  are  also  listed  in  Table  1 1 . 

The  above  steps  show  a  simplified  version  of  our  algorithm.  Some  errors  can  exist,  as  in  general  the  simulated 
phase  distributions  were  less  than  the  actual  values  due  to  the  application  of  the  32  x  32  high-pass  filter. 
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Additional  background  phase  in  images  unremoved  by  the  high-pass  filter  can  significantly  reduce  the 
quantified  magnetic  moment  of  each  straw.  However,  increasing  the  high-pass  filter  size  will  also  reduce  the 
quantified  magnetic  moment.  This  is  the  reason  why  this  algorithm  needs  to  be  applied  to  a  shorter  echo  time 
first  and  then  to  a  long  echo  time  from  a  multiple  echo  gradient  echo  sequence.  Results  in  Table  1 1  for  straws 
#3  and  #4  were  from  TE  =  15.24  ms  in  step  1,  while  for  straw  #2,  it  was  from  TE  =  10.46  ms  with  a 
multiplication  of  3  (rather  than  2)  in  step  4.  For  straw  #1,  TE  =  8.06  ms  should  be  used  in  step  1  and  then  a 
multiplication  factor  of  3  or  4  can  be  used.  However,  in  this  case  the  uncertainty  will  be  amplified  due  to  the 
multiplication  factor. 

8)  Susceptibility  quantification  of  calcium 

We  purchased  CaCI2  in  powder  forms  from  sigma.com  and  prepared  them  to  be  4  different  solutions.  Their 
concentrations  and  quantified  susceptibility  values  using  the  high-pass  filter  and  CISSCO  are  listed  in  Table 
13.  In  order  to  achieve  desired  susceptibility  values,  we  have  to  prepare  very  high  concentrations  of  calcium. 
The  quantified  susceptibility  values  (except  from  the  first  straw)  are  self-consistent  among  echo  times,  with 
differences  within  5%  (considered  as  uncertainty).  The  quantified  susceptibility  values  from  the  first  straw  vary 
by  roughly  10%,  due  to  the  high-pass  filter  effect  on  large  magnetic  moment  values.  Note  that  the  negative 
signs  of  susceptibility  values  indicate  that  calcium  solutions  are  diamagnetic  to  water. 


Straw  No. 

Concentration  (mg/ml_) 

TE  (ms) 

AX  (ppm) 

1 

503.9 

8.07 

-1.79 

10.46 

-1.60 

12.85 

-1.47 

2 

252.0 

8.07 

-1.13 

10.46 

-1.11 

12.85 

-1.06 

3 

126.0 

8.07 

-0.63 

10.46 

-0.60 

12.85 

-0.60 

4 

63.0 

8.07 

-0.35 

10.46 

-0.33 

12.85 

-0.32 

Table  13:  Susceptibility  quantifications  of  CaCI2  solutions  using  the  high-pass  filter  and  CISSCO.  The  first 
column  lists  the  straw  number.  The  second  column  lists  the  concentration  inside  each  straw.  The  third  column 
lists  the  echo  time  of  the  images  that  we  quantified  the  susceptibility.  The  fourth  column  lists  the  measured 
susceptibility  corresponding  to  each  echo  time. 

To  compare  results  with  those  from  the  high-pass  filter,  we  have  also  applied  SHARP  to  remove  the 
background  phase  and  results  are  listed  in  Table  14.  The  results  shown  in  Table  14  are  consistent  for  each 
straw  at  different  echo  times.  The  averaged  susceptibility  values  are  calculated  from  the  three  shorter  echo 
times,  as  at  TE  =  29.58  ms,  the  quantified  susceptibility  values  are  normally  smaller  due  to  the  use  of  SHARP. 
We  will  investigate  the  reason  in  the  future.  The  quantified  susceptibility  values  from  both  high-pass  filter  and 
SHARP  are  plotted  in  Fig.  14.  The  susceptibility  values  from  high-pass  filter  plotted  in  Fig.  14  are  chosen  from 
TE  =  8.07  ms,  which  has  the  minimal  effect  among  the  three  echo  times  listed  in  Table  13. 


Straw  No. 

Ax  (TE=8.07ms) 

Ax  (TE=1 0.46ms) 

Ax  (TE=1 2.85ms) 

Ax  (TE=29.58ms) 

Ax  (avg.) 

1 

-1.93 

-1.90 

-1.92 

-1.91 

-1.92 

2 

-1.16 

-1.15 

-1.16 

-1.13 

-1.16 

3 

-0.63 

-0.62 

-0.61 

-0.61 

-0.62 

4 

-0.35 

-0.34 

-0.34 

-0.32 

-0.34 

Table  14:  Susceptibility  quantifications  of  calcium  solutions  using  SHARP  and  CISSCO.  The  first  column  lists 
the  straw  number.  The  second  to  the  fifth  column  list  the  susceptibility  quantified  at  the  given  echo  time.  The 
last  column  lists  the  averaged  susceptibility  values  from  the  three  shorter  echo  times. 


22 


We  have  further  checked  the  phase  values  inside  straws  for  the  calcium  solutions,  against  the  phase  values 
calculated  from  quantified  susceptibility  values  listed  in  Tables  13  and  14.  Table  15  shows  that  the  phase 
values  after  using  SHARP  to  remove  the  background  phase  are  closer  to  and  thus  more  consistent  with  phase 
values  calculated  from  quantified  susceptibilities.  For  easier  comparisons,  some  phase  values  are  unwrapped. 
As  the  inverse  of  SNR  is  one  standard  deviation  of  the  noise  in  phase  images,  the  last  column  of  Table  15 
allows  us  to  see  whether  the  differences  between  measured  and  expected  phase  values  are  due  to  the 
presence  of  noise.  It  is  clear  that  the  differences  are  not  only  due  to  noise  but  also  due  to  possibly  the  method 
of  SHARP  itself.  In  the  human  brain  model  section  later,  we  will  show  some  evaluations  of  the  SHARP  method. 
We  will  continue  to  investigate  and  evaluate  the  outcome  of  SHARP. 


Straw 

Echo  Time 
(ms) 

Perpendicular 

Parallel 

Ph3S6theory 

Ph3S6actual 

Ph3S6theorv 

Ph3S6actual 

1/SNR 

SHARP 

1 

8.07 

2.01 

1.59 

-4.02 

-4.14 

0.024 

10.46 

2.56 

2.03 

-5.11 

-5.37 

0.024 

12.85 

3.18 

2.48 

-6.35 

-5.85 

0.025 

2 

8.07 

1.21 

0.72 

-2.41 

-2.41 

0.042 

10.46 

1.55 

1.00 

-3.10 

-3.01 

0.042 

12.85 

1.91 

1.20 

-3.83 

-4.02 

0.042 

3 

8.07 

0.66 

0.28 

-1.32 

-1.28 

0.025 

10.46 

0.83 

0.39 

-1.66 

-1.65 

0.024 

12.85 

1.01 

0.53 

-2.03 

-2.02 

0.025 

4 

8.07 

0.37 

0.06 

-0.73 

-0.65 

0.040 

10.46 

0.46 

0.15 

-0.91 

-0.84 

0.040 

12.85 

0.56 

0.29 

-1.11 

-1.00 

0.042 

High-Pass 

Filter 

1 

8.07 

1.86 

1.14 

2.55 

1.94 

0.024 

10.46 

2.16 

1.33 

1.97 

0.66 

0.024 

12.85 

2.43 

1.35 

1.42 

-0.24 

0.025 

2 

8.07 

1.18 

0.55 

-2.36 

-2.24 

0.042 

10.46 

1.50 

0.73 

-2.99 

N/A 

0.042 

12.85 

1.75 

0.70 

2.77 

2.37 

0.042 

3 

8.07 

0.66 

0.27 

-1.31 

-0.99 

0.025 

10.46 

0.81 

0.37 

-1.63 

-1.16 

0.024 

12.85 

0.99 

0.37 

-1.98 

-1.26 

0.025 

4 

8.07 

0.36 

0.05 

-0.73 

-0.5 

0.040 

10.46 

0.45 

0.12 

-0.90 

-0.59 

0.040 

12.85 

0.53 

0.10 

-1.07 

-0.69 

0.042 

Table  15:  Comparisons  between  measured  and  expected  phase  values  inside  each  straw  for  calcium  solutions 
at  both  perpendicular  and  parallel  orientations  at  three  different  echo  times.  The  first  column  lists  the  method  of 
choice  (SHARP  or  high-pass  filter)  to  remove  the  background  field.  The  second  column  lists  the  straw  number. 
The  third  column  lists  the  echo  time  of  the  images  for  quantifications.  The  fourth  and  sixth  column  list  the 
expected  phase  values  calculated  from  quantified  susceptibility  values  listed  in  Tables  13  and  14.  Note  that 
susceptibility  values  from  Tables  13  and  14  are  used  for  calculating  phase  values  based  on  the  use  of  the 
high-pass  filter  or  SHARP,  respectively.  The  fifth  and  seventh  column  list  the  measured  phase  values.  The  last 
column  lists  the  noise  of  one  standard  deviation  in  phase  images  for  the  parallel  orientation.  Signal-to-noise 
ratios  were  measured  for  each  straw  at  different  echo  times.  N/A  in  the  table  means  no  phase  value  was 
measured  due  to  reasons  such  as  obvious  artifacts. 
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Figure  14:  Quantified  susceptibility  values  based  on  CISSCO  with  either  the  use  of  the  high-pass  filter  or 
SHARP  versus  concentrations  of  CaCI2.  Linear  relations  between  the  susceptibility  and  concentration  are 
given  in  the  chart.  This  figure  clearly  shows  that  SHARP  leads  to  better  results  than  the  high-pass  filter.  Note 
that  the  slopes  of  the  fitted  lines  are  negative  as  CaCI2  is  diamagnetic  to  water. 

9)  Susceptibility  quantification  of  ferritin 

Similar  to  the  calcium  case,  we  have  quantified  the  susceptibility  of  4  different  ferritin  concentrations.  Table  16 
lists  the  results  from  the  high-pass  filter.  Again,  quantified  susceptibility  values  from  different  echo  times  are 
only  consistent  at  low  concentrations.  The  high-pass  filter  leads  to  inconsistent  values  at  high  concentrations. 


Straw  No. 

Concentration  (mg/mL) 

TE  (ms) 

Ay  (ppm) 

1 

1.87 

8.07 

2.25 

10.46 

2.08 

12.85 

1.97 

2 

0.93 

8.07 

1.25 

10.46 

1.22 

12.85 

1.12 

3 

0.47 

8.07 

0.63 

10.46 

0.63 

12.85 

0.60 

4 

0.23 

8.07 

0.32 

10.46 

0.33 

12.85 

0.32 

Table  16:  Quantified  susceptibility  values  for  ferritin  solutions  with  the  use  of  the  high-pass  filter.  The  second 
column  lists  the  iron  concentrations  rather  than  the  ferritin  concentrations.  The  meanings  of  other  columns 
have  been  explained  in  Table  13. 

To  compare  quantified  susceptibility  values  from  CISSCO,  between  the  uses  of  the  high-pass  filter  and  SHARP 
to  remove  the  background  field,  we  obtain  results  from  SHARP  and  show  them  in  Table  17.  When  ferritin 


24 


concentrations  are  relatively  low,  we  see  from  Tables  16  and  17  that  quantified  susceptibility  values  agree  with 
each  other,  regardless  of  whether  the  high-pass  filter  or  SHARP  is  applied  to  original  phase  images.  At  a  high 
concentration,  SHARP  (from  Table  17)  shows  more  consistent  results  than  the  high-pass  filter  (from  Table  16). 
As  the  results  shown  in  Table  17  are  very  consistent  for  each  straw  at  the  first  three  echo  times,  the  averaged 
susceptibility  values  are  calculated  only  from  the  three  shorter  echo  times.  At  TE  =  29.58  ms,  the  quantified 
susceptibility  values  are  slightly  smaller  due  to  the  use  of  SHARP.  The  quantified  susceptibility  values  from 
both  high-pass  filter  and  SHARP  are  plotted  in  Fig.  15.  The  fit  from  SHARP  shows  a  good  linear  relation 
between  susceptibility  and  sample  concentration.  The  slope  of  the  fit  is  1.38  ±  0.01  ppm*ml/mg,  which  is  more 
than  20%  higher  than  our  previous  results  [Zheng  et  al.,  2013].  This  difference  is  likely  due  to  the  use  of 
different  ferritin  samples  and  it  is  for  this  conjecture  we  are  re-doing  experiments  on  ferritin  solutions.  In 
addition,  our  uncertainty  0.01  ppm*ml/mg  is  likely  underestimated,  as  error  due  to  each  susceptibility 
measurement  from  CISSO  and  SHARP  was  not  included  in  the  uncertainty  estimations.  The  susceptibility 
values  from  high-pass  filter  plotted  in  Fig.  15  are  again  chosen  from  TE  =  8.07  ms,  which  has  the  minimal 
effect  among  the  three  echo  times  listed  in  Table  16. 


Straw  No. 

Ax  (TE=8.07ms) 

Ax  (TE=1 0.46ms) 

Ax  (TE=12.85ms) 

AX  (TE=29.58ms) 

Ax  (avg.) 

1 

2.59 

2.59 

2.58 

2.53 

2.59 

2 

1.30 

1.32 

1.31 

1.30 

1.31 

3 

0.65 

0.65 

0.63 

0.63 

0.64 

4 

0.33 

0.33 

0.32 

0.31 

0.33 

Table  17:  Quantified  susceptibilities  of  ferritin  solutions  using  SHARP  and  CISSCO.  The  meaning  of  each 
column  has  been  explained  in  Table  14. 


Figure  15:  Quantified  susceptibility  values  based  on  CISSCO  with  either  the  use  of  the  high-pass  filter  or 
SHARP  versus  ferritin  concentrations.  Linear  relations  between  the  susceptibility  and  concentration  are  given 
in  the  chart.  This  figure  clearly  shows  that  SHARP  leads  to  more  consistent  results  than  the  high-pass  filter. 

We  have  again  checked  the  phase  values  inside  straws  for  the  ferritin  solutions,  against  the  phase  values 
calculated  from  quantified  susceptibility  values  listed  in  Tables  16  and  17.  The  results  are  listed  in  Table  18. 
For  easier  comparisons,  some  phase  values  have  been  unwrapped  due  to  the  subtle  details  in  phase 
unwrapping  software.  We  also  list  the  inverse  of  SNR  inside  each  straw  in  Table  18  for  the  parallel  orientation, 
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as  1/SNR  is  one  standard  deviation  of  the  noise  in  phase  images.  However,  in  most  cases,  the  measured 
phase  values  differ  from  expected  phase  values  by  more  than  two  standard  deviations.  This  suggests  that 
other  scientific  reasons  rather  than  the  Gaussian  noise  are  needed  to  explain  those  differences. 


Straw 

Echo 

Time  (ms) 

Perpendicular 

Parallel 

Ph3S6theorv 

Ph3S6actual 

Phl3S6theory 

Ph3S6actual 

1/SNR 

SHARP 

1 

8.07 

-2.69 

-2.27 

5.39 

5.43 

0.020 

10.46 

-3.49 

-3.11 

6.99 

7.46 

0.021 

12.85 

-4.26 

-3.77 

8.53 

8.85 

0.024 

29.58 

-9.65 

-8.44 

19.31 

20.82 

0.050 

2 

8.07 

-1.35 

-1.54 

2.70 

2.91 

0.026 

10.46 

-1.77 

-1.95 

3.55 

3.82 

0.028 

12.85 

-2.17 

-2.44 

4.33 

4.69 

0.029 

29.58 

-4.97 

-5.69 

9.94 

10.74 

0.042 

3 

8.07 

-0.67 

-0.71 

1.35 

1.48 

0.024 

10.46 

-0.88 

-0.91 

1.76 

1.96 

0.024 

12.85 

-1.05 

-1.14 

2.10 

2.38 

0.024 

29.58 

-2.40 

-2.83 

4.80 

5.55 

0.029 

4 

8.07 

-0.34 

-0.33 

0.68 

0.80 

0.017 

10.46 

-0.44 

-0.43 

0.89 

1.08 

0.017 

12.85 

-0.54 

-0.55 

1.08 

1.26 

0.017 

29.58 

-1.20 

-1.37 

2.40 

2.86 

0.019 

High-Pass 

Filter 

1 

8.07 

-2.34 

-2.01 

-1.60 

-0.41 

0.020 

10.46 

-2.80 

-2.54 

-0.68 

0.99 

0.021 

12.85 

3.02 

N/A 

2.23 

2.27 

0.024 

2 

8.07 

-1.30 

-1.53 

2.60 

2.98 

0.026 

10.46 

-1.65 

-1.95 

-2.99 

-2.22 

0.028 

12.85 

-1.85 

-2.47 

-2.57 

-1.24 

0.029 

3 

8.07 

-0.66 

-0.71 

1.32 

1.41 

0.024 

10.46 

-0.84 

-0.85 

1.69 

1.80 

0.024 

12.85 

-0.99 

-1.01 

1.99 

2.08 

0.024 

4 

8.07 

-0.34 

-0.43 

0.68 

0.83 

0.017 

10.46 

-0.44 

-0.55 

0.88 

0.96 

0.017 

12.85 

-0.53 

-0.69 

1.07 

1.02 

0.017 

Table  18:  Comparisons  between  measured  and  expected  phase  values  inside  each  straw  for  ferritin  solutions 
at  both  perpendicular  and  parallel  orientations  at  different  echo  times.  The  meaning  of  each  column  is 
described  in  the  caption  of  Table  15,  except  that  the  expected  phase  values  (Phasetheory)  are  calculated  from 
results  in  Tables  16  or  17,  depending  on  whether  SHARP  or  the  high-pass  filter  was  used  to  remove  the 
background  phase. 

10)  Status  of  SQUID-based  magnetometer  for  measuring  the  magnetic  moment  of  each  sample 

Due  to  the  national  liquid  helium  crisis,  we  have  not  been  able  to  measure  our  samples  with  a  SQUID-based 
magnetometer  in  our  Physics  Department.  In  responding  to  the  helium  crisis,  the  Physics  Department  had 
already  moved  all  helium  consumed  equipment  to  the  basement  and  reinstalled  them  between  April  and  July 
2013.  However,  after  the  reinstallations,  the  Physics  Department  had  identified  some  problems  in  the 
basement.  Thus  they  are  now  in  the  process  of  moving  all  the  helium  consumed  equipment  to  the  third  floor  of 
their  building.  We  are  hoping  that  they  will  be  able  to  reinstall  all  the  equipment  by  the  end  of  this  year  and  then 
we  will  be  able  to  measure  the  magnetic  moment  of  each  sample  with  the  SQUID-based  magnetometer. 
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III.  Human  brain  modeling 


Our  intended  human  brain  model  (see  Fig.  16)  includes  the  red  nucleus,  substantia  nigra,  crus  cerebri, 
thalamus,  caudate  nucleus,  putamen,  globus  pallidus,  vessels,  grey  matter,  white  matter  in  the  cerebral  cortex, 
and  the  cerebellum  along  with  cerebrospinal  fluid  (CSF).  This  model  is  a  generalized  model  of  the  structures. 
The  input  susceptibility  value  of  each  tissue  is  given  in  Table  19  and  other  tissue  parameters  are  shown  in  Fig. 
17.  Relative  positions  of  the  structures  will  vary  for  different  individuals  and  this  inter-individual  variability  is 
expected.  The  knowledge  of  the  general  shape  and  positions  of  the  structures  in  the  brain  model  helps  to  study 
a  general  phase  behavior  of  the  structures  relative  to  each  other  in  the  same  3D  space. 


Figure  16:  Transverse,  Coronal,  and  Sagittal  views  of  the  3D  brain  model.  The  structures  are  differentiated  by 
their  susceptibility  values  which  are  listed  in  Table  19.  CN  =  Caudate  nucleus,  SN  =  Subtantia  nigra,  RN  =  Red 
nucleus,  GP  =  Globus  pallidus,  CC  =  Crus  Cerebri,  PUT  =  Putamen,  TH  =  Thalamus,  WM  =  White  matter,  GM 
=  Grey  matter,  CSF  =  Cerebrospinal  fluid. 


Structure 

X  value  (in  ppm) 

Structure 

X  value  (in  ppm) 

White  matter 

0 

Veins 

0.45 

Grey  matter 

0.02 

Red  Nucleus 

0.13 

Globus  Pallidus 

0.18 

Substantia  Nigra 

0.16 

Putamen 

0.09 

Thalamus 

0.01 

Caudate  Nucleus 

0.06 

Crus  Cerebri 

-0.03 

Cerebrospinal  Fluid 

-0.014 

Table  19:  Susceptibility  values  used  in  our  human  brain  model. 
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Figure  17:  Spin  density  (p0),  Ti  and  T2*  properties  are  added  to  each  tissue  in  the  brain  model  in  order  to 
generate  magnitude  images  at  main  magnetic  field,  B0,  of  3  Tesla  for  different  echo  times  (TE  =  4  ms,  7  ms,  10 
ms,  13  ms,  16  ms,  19  ms,  22  ms,  and  25  ms)  using  Eq.  (2)  where  flip  angle  of  15°  is  used.  The  reasonable 
values  for  p0,  Ti  and  T2*  are  selected  from  real  data  processing  (for  spin  density  and  T2*  values  of  basal 
ganglia  structures)  and  literature  at  TR  =  35  ms. 


1)  Simulations  of  magnitude  images  and  phase  images 


We  first  simulate  the  brain  with  a  matrix  size  of  512  x  512  x  512  grid  points,  which  corresponds  to  an  image 
resolution  of  0.5  X  0.5  X  0.5mm3.  The  biological  tissues  used  in  the  model  have  different  relaxation  times  T-i 
and  T2‘,  given  in  Fig.  17,  and  susceptibility  values  listed  in  Table  19. 


The  magnitude  signal-intensity  from  an  rf-spoiled  short-TR  gradient-echo  sequence  is  given  by 


S  =  p^sinQ  exp  (-  X  1  -  exp (-  ^}]  /[I  -  exp (- 


(2) 


where  TR  is  the  repetition  time  and  0  is  the  flip  angle. 

2)  Simulations  of  partial  volume  effects  for  tissues  including  small  and  large  vessels 

The  phase  simulations  of  the  brain  model  can  be  compared  with  the  in  vivo  human  phase  images  by 
introducing  confounding  factors,  such  as  partial  volume  effects  (due  to  discretization  of  the  MR  signal)  and 
white  Gaussian  noise  (incorporated  through  real  and  imaginary  part  of  signals). 

In  order  to  simulate  the  effect  due  to  partial  volume  and  Gibbs  ringing,  we  performed  a  process  analogous  to 
the  MRI  image  acquisition.  We  start  by  simulating  magnitude  and  phase  images,  with  a  matrix  size  of  2048  x 
2048  x  2048  grid  points.  Then,  lower  resolution  images  are  obtained  by  taking  the  central  128  x  128  x  128 
points  of  the  Fourier  transform  of  the  high  density  matrix,  followed  by  an  inverse  Fourier  transform.  This  new 
data  comprises  of  experimental  artifacts  such  as  Gibbs  ringing  and  partial  volume  effects ,  making  the  phase 
behavior  much  more  realistic. 
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Figure  18  shows  the  changes  in  susceptibility  values  due  to  the  partial  volume  effect,  before  and  after 
applications  of  post-processing  techniques  such  as  homodyne  high-pass  filter  (with  different  sizes)  and 
SHARP.  The  susceptibility  values  were  quantified  using  our  developed  SWIM  (susceptibility  weighted  imaging 
mapping)  program  [Haacke  et  al.,  2010].  The  diameter  of  a  septal  vein  (2  pixels)  is  much  smaller  than  that  of 
the  straight  sinus  (8  pixels).  Hence,  the  phase  integration  due  to  the  partial  volume  effect  clearly  affects  the 
susceptibility  quantification  of  septal  vein  more  than  that  of  the  straight  sinus. 
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Figure  18:  Mean  susceptibility  value  measurements  of  (a)  the  straight  sinus  and  (b)  septal  veins,  from  the 
susceptibility  maps,  generated  from  filtered  and  unfiltered  phase,  before  and  after  simulating  the  partial  volume 
effects  (and  mimicking  the  discrete  sampling  of  the  MR  signal).  The  susceptibility  value  of  the  septal  vein  is 
severely  affected  by  the  partial  volume  effect.  The  diameter  of  the  straight  sinus  is  around  8  pixels  and  the 
diameter  of  the  septal  vein  is  2  pixels  in  the  model.  The  initial  susceptibility  value  inside  these  veins  was  set  to, 
Ax  =  0.45  ppm. 

3)  Estimating  the  ideal  echo  time  for  SWIM 

One  of  the  critical  features  in  assessing  iron  with  quantitative  susceptibility  mapping  (QSM)  is  the  choice  of 
echo  time.  This  is  particularly  important  for  traumatic  brain  injury  where  high  iron  content  can  exist  in 
microbleeds. 

As  seen  in  Fig.  19,  the  error  in  the  susceptibility  maps  decreases  as  we  increase  the  echo  time.  The  reduction 
in  error  seems  significant  between  TE  =  4  ms  and  TE  =  13  ms.  It  should  be  noted  that  when  the  echo  time 
equals  to  the  T2*  value  of  the  structure,  it  produces  high  SNR  in  phase  images.  The  standard  deviation  in 
phase  images  is  equal  to  the  inverse  of  SNR  in  magnitude  images.  This  may  be  the  reason  why  we  see  the 
lowest  standard  deviation  at  around  TE  =  20  ms  (T2*(giobuspaiiidus)  =  19  ms  in  the  model).  The  susceptibility  values 
generated  from  echo  times  longer  than  16  ms  do  not  show  much  reduction  in  errors.  This  fact  suggests  that 
using  echo  time  of  13  ms  will  be  a  good  practice  for  susceptibility  quantifications  with  SWIM,  as  such  a  choice 
allows  us  to  shorten  the  acquisition  time. 
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Error  in  /  quantification  in  homogeneous  white  matter 


a) 


Mean  /  and  standard  deviation  quantflcation  inside  Globus  pallidus 


b) 


Figure  19:  The  relationship  between  echo  time  and  error  in  susceptibility  quantifications  using  SWIM.  A 
Gaussian  noise  (mean/standard  deviation  =  10/1)  was  added  to  the  real  and  imaginary  parts  of  the  simulated 
data,  (a)  The  measurements  are  performed  in  a  homogeneous  white  matter  region  where  the  expected  mean 
value  is  ideally  zero.  A  small  homogeneous  white  matter  region  of  50  pixels  in  size  was  chosen  to  calculate  the 
mean  susceptibility  value,  (b)  The  measurements  of  mean  %  value  and  the  standard  deviation  are  calculated 
inside  the  Globus  pallidus  (x  (giobus  pamdus) =  0.18  ppm).  A  total  of  100  pixels  inside  the  globus  pallidi,  avoiding  the 
edges,  were  used  for  calculating  the  mean  susceptibility  and  its  associated  uncertainty. 

4)  Human  brain  model  simulations  with  more  grid  points  and  results 

We  have  further  properly  improved  and  simulated  human  brain  magnitude  and  phase  images  by  including  the 
partial  volume  effects  due  to  discrete  voxels  containing  MR  signals,  with  a  final  matrix  size  of  2563.  This  matrix 
size  mimics  the  typical  1  mm  isotropic  images  for  better  depiction  of  tissues.  This  simulation  procedure 
required  a  very  long  computer  time,  as  we  first  generated  a  40963  complex  matrix,  Fourier  transformed  the 
40963  complex  matrix,  took  the  central  2563  points  of  the  matrix,  and  performed  an  inverse  Fourier  transform  of 
the  2563  matrix  to  obtain  images  in  the  spatial  domain.  Such  a  ratio  (or  higher)  between  the  two  matrices  is 
needed  for  properly  generating  the  partial  volume  effect  for  relatively  small  objects.  The  Gaussian  noise  can  be 
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added  to  the  2563  complex  images.  This  new  dataset  comprises  of  actual  artifacts  such  as  Gibbs  ringing  and 
partial  volume  effects,  making  the  phase  behavior  much  more  realistic.  Examples  are  shown  in  Fig.  20.  We 
can  see  clear  differences  of  small  structures  between  the  high  and  low  density  matrices. 


Figure  20:  Comparisons  between  the  simulated  phase  images  without  partial  volume  effects  (a)  and  (c)  and 
with  the  partial  volume  effects  (b)  and  (d).  Imaging  parameters  used  for  (a)  and  (b)  are  B0  =  3  T  and  TE=  15  ms 
and  for  (c)  and  (d)  are  B0  =  3  T  and  TE  =  29  ms.  Partial  volume  effects  are  generated  by  cropping  the  central 
256  x  256  x  256  k-space  elements  from  a  high  grid  matrix  of  4096  x  4096  x  4096.  K-space  cropping  also 
introduces  Gibbs  ringing  in  the  image. 

Addition  of  the  partial  volume  effects  in  the  simulations  is  a  key  factor  in  understanding  the  actual  phase  data. 
When  the  current  SWIM  technique  was  applied,  the  mean  susceptibility  inside  a  relatively  small  vessel,  such 
as  the  right  septal  vein  with  a  diameter  of  2  pixels  (white  arrows  in  Fig.  21),  shows  significant  reduction  of  the 
susceptibility  value  after  introducing  the  partial  volume  effect.  The  quantified  susceptibility  inside  the  right 
septal  vein  quantified  by  SWIM  based  on  original  images  without  the  partial  volume  effect  is  Ax  =  0.38  ±  0.02 
ppm,  which  is  about  16%  lower  than  the  actual  value  of  Ax  =  0.45  ppm.  On  the  other  hand,  the  mean 
susceptibility  value  of  the  same  vein  quantified  by  SWIM  from  images  with  the  partial  volume  effect  is  Ax  = 
0.17  ±  0.03  ppm,  with  a  much  worse  error  of  62%.  These  values  also  agree  with  results  shown  in  Fig.  18. 
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Figure  21:  Maximal  Intensity  Projection  (MIP)  images  of  the  susceptibility  maps,  (a)  MIP  (over  42  slices)  of  the 
susceptibility  maps  generated  from  the  original  phase  images,  and  (b)  MIP  (over  the  same  FOV  as  in  (a))  of 
the  susceptibility  maps  generated  from  phase  images  with  the  partial  volume  effect.  The  susceptibility  values 
inside  veins  are  fairly  a  constant  for  (a)  irrespective  of  the  size  of  a  vessel,  whereas  in  (b)  the  susceptibility 
values  of  veins  show  variations,  depending  on  the  vessel  size.  The  white  arrow  identifies  the  right  septal  vein, 
where  the  reduction  can  be  clearly  seen  in  (b). 

After  the  above  simulations,  Gaussian  noise  is  added  to  the  real  part  and  imaginary  part  of  the  2563  complex 
images.  An  example  image  with  an  SNR  of  20:1  is  shown  in  Fig.  22.  Constructing  such  a  numerical  phantom 
through  forward  simulations  with  known  input  susceptibility  values  is  essential  for  our  tests  of  methods  in  the 
future.  Below  we  present  some  examples. 


Figure  22:  An  example  of  the  phase  images  including  the  partial  volume  effects,  Gibbs  ringing  effects,  and 
Gaussian  noise. 

5)  Accuracy  in  background  phase  removal  methods  with  SWIM:  SHARP,  homodyne  high-pass  filtering, 
and  variable  high-pass  filtering  tested  on  our  human  brain  model 

In  this  section,  the  errors  in  estimated  susceptibility  values  using  different  background  phase  removal  methods 
are  evaluated  on  our  simulated  3D  brain  model.  Particularly,  three  phase  processing  methods  are  evaluated: 
SHARP,  homodyne  high-pass  filtering  (HP)  and  variable  high-pass  filtering  (VHP).  Homodyne  high-pass 
filtering  is  the  most  traditional  phase  processing  method  and  it  and  SHARP  have  been  also  used  with  CISSCO 
in  earlier  sections.  It  is  known  that  high-pass  filtering  process  will  lead  to  reduction  of  phase  values,  especially 
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to  relatively  large  structures.  SHARP  is  based  on  the  spherical  mean  value  property  of  the  background  field, 
but  the  accuracy  of  SHARP  depends  on  the  choices  of  the  kernel  size  (i.e.,  volume  of  a  sphere  inside  which 
phase  values  are  added)  and  the  regularization  threshold.  The  variable  high-pass  filtering  is  also  based  on 
spherical  mean  value  theorem.  It  uses  variable  filter  sizes  to  reduce  the  signal  loss  to  structures  in  central  part 
of  the  brain. 

a)  Theory 

SHARP 

Assuming  that  the  total  field  variation  &&  can  be  considered  as  a  combination  of  the  local  field  variation  and 
the  background  field  variation 

AB  =  ARi-k-AB^  [3], 

In  different  phase/field  processing  methods,  different  properties  of  the  background  field  are  utilized.  In  SHARP 
it  is  assumed  that  the  background  field  variation,  AEu  has  the  spherical  mean  value  property  (SMV)  inside  the 
brain  region  (the  value  of  certain  pixel,  e.g.,  A,  of  is  the  mean  field  value  for  all  pixels  in  a  sphere, 
centered  at  A). 

The  first  step  of  SHARP  is  similar  to  a  high-pass  filtering 


AB  —  AB  *  i  ■  -  AB ;  *  f)  "  C.AB jj  -  AB &  *  ■  AB ;  —  AB ;  *  s  [4]i 

where  s  is  a  normalized  sphere  in  image  domain  (i.e.,  a  numerical  sphere  which  is  1/num  for  pixels  inside  the 
sphere  and  0  outside,  and  “num”  is  the  number  of  pixels  inside  the  sphere.  Thus,  this  is  essentially  an 
averaging  filter.  *  denotes  for  the  convolution  process). 

Let’s  denote  the  result  of  Eq.  4  as  AB' .  Since  the  convolution  does  not  give  correct  value  near  the  edge,  an 
eroded  mask,  M  has  to  be  applied.  Equation  4  with  the  application  of  the  mask  can  be  written  as 

AU'-ADi*  (VS  -  [5], 

It  is  possible  to  solve  Eq.  5  for  AS; 

AS ■  =  AS'  *  (cf  -  S’)"1  [6] 

where  e?  -  -T1  is  the  inverse  of  the  kernel  (5  -  S).  This  inverse  process  has  to  be  regularized  typically 
through  a  Fourier  domain  truncation  (TSVD)  by  setting  a  threshold,  th.  i.e,.  the  inverse  of  the  Fourier  transform 
of  (5  -  S)  is  set  to  0  when  the  absolute  value  of  the  Fourier  transform  of  (5  -  S)  is  below  th.  Obviously,  the 
accuracy  of  SHARP  is  dependent  on  the  size  of  the  sphere  that  is  used  in  Eq.  4,  and  the  threshold  th  is  the 
deconvolution  process  in  Eq.  6. 

Homodyne  high-pass  filtering 

In  Homodyne  high-pass  filtering,  the  background  field  component  in  Eq.  3  is  assumed  to  have  low  spatial 
frequency  and  hence  can  be  removed  through  a  high-pass  filtering,  as  shown  below 

Mct& ■  FT~l  [F7[Ma&  ■  ■  Harm)  [7] 

=  *•*(<*  [8], 

The  size  of  the  Hann  window  in  k-space  in  Eq.  7  is  usually  selected  as  64  x  64  in  order  to  obtain  the  low 
spatial-frequency  components  and  to  avoid  too  much  signal  loss  to  the  local  field  variation. 
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b)  Theory:  variable  high-pass  filtering 

In  variable  high-pass  filtering,  the  background  field  component  will  be  estimated  by  applying  a  spherical  mean 
filter  to  the  original  phase  images 

AB  *s  =  AS;,  -F  AS;  [9], 

and  the  high-pass  filtered  image  is  calculated  as 

AEfrf  ■  AB  —  AB  fi  j  j-jqj 

This  is  the  same  as  Eq.  4.  However,  when  the  size  of  s  is  large  enough,  the  average  value  of  the  local  field  is 
close  to  0  (A|CJ  ).  Thus, 

ABfri?  v  AB  “  AB ■  AB •  m  -m 


However,  the  size  of  the  spherical  kernel  s  is  limited  by  the  dimension  of  the  image  matrix.  In  addition,  since 
the  spherical  mean  value  property  cannot  be  applied  across  the  edge  of  the  brain,  we  gradually  reduce  the 
size  of  the  spherical  filter  from  the  central  part  of  the  brain  to  the  periphery.  A  large  spherical  window  is  used 
for  the  central  region  (and  is  applied  in  image  domain  and  can  be  considered  as  a  very  mild  high-pass  filter), 
and  smaller  spherical  windows  (stronger  high-pass  filter)  are  used  for  regions  close  to  the  edge. 

c)  Method 

Phase  simulations  and  processing 


The  phase  images  of  the  3D  brain  model  were  calculated  using  the  fast-forward  field  calculation  shown  in 
Cheng  et  al.  [2009b]  with  a  magnetic  field  strength  of  3  T  at  an  echo  time  of  10  ms.  The  susceptibility  of  the 
sinuses  was  set  to  be  9  ppm,  in  order  to  simulate  the  background  field  variations  induced  by  the  air-tissue 
interfaces.  Phase  images  of  the  brain  structures  $;.5tn;and  of  the  sinuses  were  calculated 

independently.  The  sum  of  these  two  components  is  referred  to  as  the  phase  images  of  the  brain  model 
).  An  example  image  is  shown  in  Fig.  23. 


Figure  23:  (a)  Simulated  phase  image  of  the  brain  model  and  (b)  the  reference  region  used  for  evaluating  the 
accuracy  in  the  processed  phase  images  after  the  background  phase  is  removed. 


SHARP  with  different  kernel  sizes  ranging  from  1  pixel  to  16  pixels  (i.e.,  radius  of  the  sphere  used  in  SHARP) 
and  different  thresholds  (th  ranging  from  0.005  to  0.1  with  a  step  size  0.005)  was  applied  to  $fa».  A  brain 
mask  was  used  in  SHARP  (Eq.  5).  Note  that  the  mask  has  to  be  eroded  by  different  distances  depending  on 
the  size/radius  of  the  spherical  kernel.  Thus  the  number  of  eroded  brain  pixels  was  also  measured  for  different 
kernel  sizes. 
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Homodyne  high-pass  filter  with  a  k-space  window  size  of  64  x  64  was  also  applied  to  sps-im  for  comparisons. 
Variable  high-pass  filter  was  applied  with  a  big  sphere  (radius  =  32  pixels)  for  central  regions  and  with  smaller 
spheres  for  regions  close  to  the  edge.  At  the  edge,  the  kernel  size/radius  of  the  sphere  is  1  pixel.  The  size  of 
the  spherical  kernel  for  a  particular  pixel  was  determined  as  the  distance  from  that  pixel  to  the  nearest  edge  of 
the  3D  brain. 


The  processed  phase  images,  which  correspond  to  the  local  phase  induced  by  various  brain  structures,  were 
compared  with  the  simulated  phase  of  purely  various  structures  inside  the  brain  RMSE  (defined  in  Eq. 

12)  was  calculated  in  a  reference  region  which  corresponds  to  the  eroded  mask  in  SHARP  when  the  largest 
sphere  was  used  (see  Fig.  23),  for  all  three  background  phase  removal  methods. 


RMSE  = 


[12], 


where  N  is  the  total  number  of  pixels  in  the  reference  mask. 


Susceptibility  quantification 

Susceptibility  maps  were  generated  using  a  truncated  k-space  division  approach,  with  a  threshold  value  of  0.1 
[Haacke  et  al.,  2010],  The  mean  susceptibility  values  were  measured  for  a  total  of  9  different  structures,  as 
listed  in  Table  20.  In  order  to  independently  analyze  the  errors  in  susceptibility  quantifications  induced  by 
background  phase  removal  methods,  reference  susceptibility  maps  were  generated  using  the  ideal  phase 
images  containing  only  the  local  structures,  with  the  same  k-space  threshold  value.  The  errors  in  the  reference 
susceptibility  maps  were  purely  caused  by  the  inverse  filter  in  SWIM. 

Xss:  ~  fea? 

The  relative  error  in  susceptibility  quantification  of  each  structure  is  calculated  as  Lfowl  ,  where  Kesris 
the  estimated  susceptibility  value  and  is  the  true  susceptibility  value.  The  effects  of  the  kernel  size  and 
threshold  in  SHARP  were  studied  using  the  overall  RMSE  in  the  processed  phase  and  susceptibility  maps,  as 
well  as  the  relative  error  in  estimated  susceptibility  for  each  structure.  The  overall  RMSE  in  susceptibility  maps 
were  also  calculated  in  a  similar  way  to  Eq.  12. 


d)  Results 


Methods 

True 

Susceptibility 

Original 

Phase 

SHARP  (r=8, 
th=0.02) 

SHARP 
(r=8,  th=0.05) 

HP64 

VHP 

Veins 

0.45 

0.41  ±0.04 

0.40  ±  0.04 

0.39  ±  0.04 

0.33  ±0.07 

0.40  ±  0.05 

Thalamus 

0.01 

0.01  ±0.02 

0  ±  0.02 

0  ±  0.02 

0  ±  0.02 

0  ±  0.02 

Red  Nucleus 

0.13 

0.12  ±0.01 

0.11  ±0.01 

0.11  ±0.01 

0.07  ±0.02 

0.11  ±0.01 

Substantia  Nigra 

0.16 

0.15  ±0.02 

0.14  ±0.02 

0.14  ±0.02 

0.10  ±0.04 

0.13  ±0.02 

Subthalamic 

Nucleus 

0.20 

0.19  ±0.02 

0.19  ±0.02 

0.19  ±0.02 

0.14  ±0.03 

0.18  ±0.02 

Crus  Cerebri 

-0.03 

-0.03  ±  0.02 

-0.04  ±0.02 

-0.04  ±  0.02 

-0.03  ±  0.03 

-0.04  ±  0.02 

Caudate 

0.06 

0.05  ±0.01 

0.04  ±0.01 

0.04  ±0.01 

0.02  ±  0.02 

0.04  ±0.01 

Putamen 

0.09 

0.08  ±0.01 

0.06  ±0.01 

0.04  ±0.01 

0.01  ±0.03 

0.05  ±0.01 

Globus  Pallidus 

0.18 

0.16  ±0.01 

0.14  ±0.01 

0.12  ±0.01 

0.05  ±  0.03 

0.11  ±0.01 

Table  20:  The  quantified  susceptibility  values  (mean  ±  st.  dev.  in  ppm)  using  SWIM  for  different  structures  and 
background  phase  removal  methods.  The  first  column  lists  the  9  structures.  The  second  column  lists  the  true 
susceptibility  values  in  ppm.  The  third  column  lists  results  from  the  reference  images,  without  background 
phase  added  into  images.  The  fourth  and  the  fifth  column  list  results  for  SHARP  with  radius  =  8  pixels, 
threshold  =  0.02  and  0.05,  respectively.  The  sixth  column  lists  the  results  after  the  application  of  a  64  x  64  high 
pass  filter.  The  seventh  column  lists  results  from  the  Homodyne  high-pass  filter. 
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The  RMSE  estimated  from  processed  phase  images  for  different  kernels  (radii)  and  thresholds  in  SHARP  are 
shown  in  Fig.  24a.  These  RMSE  represent  overall  error  levels  in  the  phase  images  processed  by  SHARP.  For 
a  given  kernel  size,  a  minimal  RMSE  can  be  seen  from  Fig.  24a  and  is  plotted  in  Fig.  24c.  This  minimal  RMSE 
is  decreasing  as  the  size  of  the  kernel  increases  (Fig.  24c),  but  the  information  loss  (the  number  of  eroded 
pixels)  also  increases  (Fig.  24d).  At  the  minimal  RMSE  for  a  given  kernel  size,  the  optimal  threshold  values  are 
shown  in  Fig.  24b.  As  a  comparison,  the  RMSE  for  the  HP  processed  phase  is  0.05  radian  and  the  RMSE  for 
the  VHP  processed  phase  is  0.02  radian. 
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Figure  24:  (a)  RMSE  from  phase  images  for  different  radii  of  the  spherical  kernels  and  different  threshold 
values  used  in  SHARP,  (b)  Optimal  threshold  (minimal  RMSE)  for  different  radii,  (c)  Minimal  RMSE  obtained 
from  the  optimal  threshold  shown  in  (a)  for  each  given  kernel  (radius),  (d)  Relative  information  loss  (the 
number  of  eroded  pixels  over  the  total  number  of  pixels  in  the  brain  mask)  for  different  spherical  kernels  (radii). 


The  overall  RMSE  in  the  susceptibility  maps,  generated  from  SHARP  processed  phase  images  are  shown  in 
Fig.  25a.  The  minimal  RMSE  for  each  kernel  size  is  shown  in  Fig.  25b.  These  two  figures  resemble  the  results 
shown  in  Figs.  24a  and  24c. 
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Figure  25:  (a)  Overall  RMSE  in  susceptibility  quantification  for  different  kernel  sizes  (radii)  and  different 
thresholds,  (b)  The  minimal  RMSE  for  a  given  kernel  size  (radius). 
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Generally  speaking,  larger  radii  of  spherical  kernels  lead  to  smaller  errors.  In  addition,  the  radius  of  the  sphere 
should  be  larger  than  4  pixels,  in  order  to  maintain  a  low  RMSE.  However  for  each  individual  structure  listed  in 
Table  20,  the  optimal  kernel  size  and  threshold  may  vary,  as  shown  in  Fig.  26,  which  shows  the  relative  errors 
of  susceptibility  quantifications  (in  absolute  values).  For  basal  ganglia  structures,  the  inverse  filter  itself  will 
cause  an  underestimation  around  10%-20%.  The  SHARP  and  VHP  processed  phase  images  lead  to  an 
underestimation  of  around  30%-40%.  This  suggests  that  SHARP  and  VHP  will  cause  roughly  20%  error  for 
basal  ganglia  structures.  For  other  structures  such  as  veins,  the  underestimation  is  smaller.  The  most  severe 
error  is  seen  from  thalamus.  In  this  case,  the  errors  in  the  measured  susceptibility  values  are  high  for  all  three 
phase  processing  methods.  Even  using  the  phase  images  without  the  background  field  in  susceptibility 
quantifications,  the  uncertainties  are  still  high.  This  might  suggests  that  the  problem  is  largely  caused  by  the 
inverse  filter  itself. 


From  Fig.  26,  it  can  be  seen  that  generally  the  errors  in  susceptibility  estimates  decrease  as  the  radius  of  the 
SHARP  kernel  increases  (except  for  substantia  nigra  and  crus  cerebri).  Considering  the  relatively  higher  signal 
loss  at  larger  kernel  sizes,  a  reasonable  choice  for  the  size  of  the  spherical  kernel  in  SHARP  should  be  around 
6-8  pixels.  The  threshold  value  should  be  kept  as  small  as  possible  (th  <=  0.05).  In  some  situations,  a  higher 
threshold  value  may  be  desired,  in  order  to  suppress  non-harmonic  phase  artifacts. 
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Figure  26:  Relative  errors  measured  from  susceptibilities  (in  absolute  values)  using  different  parameters  in 
SHARP  for  9  different  structures. 

In  addition,  the  estimated  susceptibilities  using  the  original  phase  images  (i.e.,  images  without  the  background 
phase)  were  plotted  against  the  estimates  using  the  processed  phase  images,  as  shown  in  Fig.  27.  Each  x- 
axis  shows  the  susceptibility  estimates  obtained  using  the  original  phase  images,  and  each  y-axis  shows  the 
estimates  obtained  through  the  three  background  phase  removal  methods.  Figures  27a  and  27b  are  obtained 
through  SHARP  with  different  regularization  parameters  (see  Table  20). 

Figure  27  demonstrates  possible  effects  of  different  background  phase  removal  methods  on  susceptibility 
quantifications.  When  using  SHARP  with  proper  thresholds,  the  estimated  susceptibilities  are  almost  the  same 
as  those  estimated  using  the  original  phase  (the  slope  is  close  to  1 ).  However,  it  is  also  noticed  that,  when  the 
threshold  value  becomes  larger,  the  correlations  between  the  estimates  become  worse.  On  the  other  hand, 
using  the  phase  images  processed  by  the  homodyne  high-pass  filter  (64  x  64  window  size),  the  susceptibilities 
are  certainly  underestimated.  Using  the  phase  images  processed  by  variable  high-pass  filter,  again  the 
estimated  susceptibilities  are  almost  the  same  as  those  obtained  from  the  original  phase  images. 
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Figure  27:  Estimated  Susceptibility  values  using  the  original  phase  vs.  different  background  phase  removal 
methods:  (a)  SHARP  (radius  =  8  pixels,  th  =  0.02),  (b)  SHARP  (radius  =  8  pixels,  th  =  0.05),  (c)  Homodyne  HP 
32  x  32,  (d)  Homodyne  HP  64  x  64,  and  (e)  Variable  HP. 

e)  Discussion 

It  can  be  concluded  that  errors  derived  from  SHARP  processed  phase  images  decrease  as  the  size  of  the 
spherical  kernel  increases  (Fig.  25).  This  is  mainly  due  to  the  deconvolution  process,  as  larger  spherical  kernel 
leads  to  less  amplification  of  any  remnant  background  field. 

Clearly  the  variable  high-pass  filtering  method  leads  to  less  underestimation  in  susceptibility  quantifications, 
compared  with  the  traditional  homodyne  high-pass  filter.  In  variable  high-pass  filter,  the  spherical  kernel  is 
used  to  generate  a  low-pass  filtered  phase  image  (the  averaged  phase  image)  which  is  close  to  the 
background  phase.  No  additional  deconvolution  is  required.  Thus  there  is  no  need  to  select  any  regularization 
parameter  as  it  is  done  in  SHARP.  As  can  be  seen  from  the  results,  variable  high-pass  filtering  leads  to  similar 
results  as  in  SHARP. 

The  SHARP  and  VHP  processed  phase  images  lead  to  small  errors  for  most  structures.  However  for  thalamus 
or  structures  with  small  susceptibility  values  (below  0.1  ppm),  even  though  some  structures  are  large  (occupy 
many  pixels),  the  errors  from  susceptibility  quantifications  are  large.  On  the  other  hand,  for  veins  and  basal 
ganglia  structures  with  relatively  larger  susceptibility  values  (above  0.1  ppm),  the  errors  seem  small.  However, 
Table  20  shows  that  the  errors  of  susceptibility  value  estimated  from  Globus  Pallidus  with  SHARP  and  VHP 
are  not  small.  This  indicates  that  some  systematic  errors  in  SHARP  and  VHP  have  not  been  identified  yet.  This 
subtlety  is  also  supported  by  results  shown  between  Table  11  and  Table  12,  as  discussed  earlier.  One 
theoretical  problem  of  SHARP  is  that  the  mean  value  theory  for  SHAPR  becomes  invalid  when  the  spherical 
kernel  crosses  the  boundary  of  an  object.  This  includes  the  situation  when  an  object  is  completely  enclosed  by 
the  volume  of  the  spherical  kernel,  as  the  mean  value  theory  was  applied  to  the  surface  of  the  spherical  kernel 
rather  than  its  volume.  Thus,  a  better  background  phase  removal  technique  is  still  useful  to  be  developed. 
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6)  Evaluating  the  errors  in  Laplacian  based  phase  processing  algorithms 

The  Laplacian  based  phase  processing  algorithms  are  a  substitute  for  the  phase  unwrapping  procedure.  As  we 
mentioned  before,  phase  unwrapping  is  needed  for  SHARP  and  quadratic  fitting  methods.  Thus,  it  becomes 
important  to  evaluate  errors  caused  by  these  algorithms.  Currently  there  are  mainly  two  ways  to  calculate  the 
Laplacian  of  phase.  One  method  is  to  use  the  phase  differences  and  discrete  Laplacian  operator,  and  the  other 
is  to  use  the  sine  and  cosine  functions  together  with  Fourier  transform.  While  the  former  method  is  essentially 
a  discrete  approximation  of  Laplacian,  the  latter  method  is  based  on  continuous  functions  of  phase.  Given  the 
fact  that  MRI  phase  images  are  discrete,  it  might  be  more  consistent  to  use  all  operators  in  discrete  forms.  In 
this  report,  the  errors  associated  with  both  phase  processing  methods  are  evaluated  using  both  numerical 
model  and  in  vivo  data. 

a)  Theory 

The  total  field  can  be  written  as  the  sum  of  the  local  field  and  background  field 
Aff(r)  =  &8b  (r)  +  AMr)  [1 3] 

In  a  homogenous  region,  the  background  field  satisfies  the  Laplace’s  equation 

[14] 

From  Eqs.  13  and  14,  the  relation  between  the  total  field  and  local  field  can  be  derived  as 
?2  [Aff(r¥]  =  Vs  [ASf  (r)]  [1  5] 

The  left  side  of  Eq.  15,  the  Laplacian  of  the  total  field  can  be  calculated  using  the  discrete  Laplacian  operator. 

+  4,  £  Wit  41  +  ft  [-|g] 

where  is  the  phase  value  of  the  pixel  at  index  (i,j,k),  and  0  =  —yABTE  for  a  right-handed  system.  The 

right  hand  side  of  Eq.  16  can  be  calculated  through  the  differences  of  two  neighboring  pixels  in  the  original 
phase  image  with  phase  wraps  (<5^).  For  example, 

&4V,1 1 "  HfrU  ■  ~  [1 7] 

In  this  study,  Eqs.  16  and  17  will  be  denoted  as  “discrete  operators”.  Alternatively,  as  proposed  by  Schofield 
and  Zhu  [2003],  the  Laplacian  of  the  unwrapped  phase  can  be  calculated  from  ^vr- 

Y s  $  -  -  3  tn<fofY s  (gqs  [  1 8] 

The  Laplacians  can  be  calculated  using  properties  of  Discrete  Fourier  Transforms 
Vs fir)  -  -  r[/{r>]}  [1 9] 


As  both  Eqs.  18  and  19  are  for  continuous  functions,  they  will  be  denoted  as  “continuous  operators”  in  this 
study.  Then  from  Eq.  15, 

Tsp  =  Vs0t,  orVS^-^>  0.  [20] 

Solving  Eq.  20  requires  boundary  condition,  which  is  usually  unknown.  However,  if  the  Laplacian  is  calculated 
using  Fourier  transform,  then  periodicity  is  imposed  on  the  boundary.  In  this  case, 
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[21] 


where  M  is  an  eroded  brain  mask,  is  the  regularized  inverse  of  1{&),  the  Fourier  transform  of  the 

discrete  Laplacian  kernel. 


'{FTtibyfi-lwkm  IffUWH  >  A 
0,  wtisn  |f  r[tC»*)]l  ^  ^ 


[22] 


where  th  is  the  truncation  threshold  in  Fourier  domain,  which  was  set  as  0.015  in  this  study.  If  the  continuous 
form  of  Laplacian  was  used,  then  the  local  phase  can  be  reconstructed  by  inverting  Eq.  19.  Regularization  is 
also  required  for  inverting  Eq.  19.  Particularly,  the  truncation  threshold  was  empirically  set  to  be  4*1  O'4.  The 
whole  process  is  illustrated  in  Fig.  28. 


Original  Phase  Phase  Laplacian 


Erosion 

Local  Phase  J 


Laplacian 

Inversion 


Figure  28:  Illustration  of  the  Laplacian  based  phase  image  processing  steps. 

b)  Method 

The  above  3D  numerical  brain  model  was  used  for  evaluating  the  Laplacian  based  algorithms.  The  background 
phase  was  added  by  including  the  air-sinuses  which  have  a  susceptibility  of  9  ppm  relative  to  the  brain  tissue. 
The  phase  images  were  simulated  with  B0  =  3  T,  TE  =  10  ms.  Aliased  phase  images  were  processed  using 
both  discrete  and  continuous  Laplacian  based  algorithms. 

The  errors  in  the  phase  processing  were  evaluated  by  comparing  the  extracted  local  phase  information,  <j>i  in 
Eq.  21,  with  the  actual  local  phase  information  produced  from  the  forward  calculations.  In  order  to  determine 
the  effects  of  phase  processing  on  susceptibility  quantifications,  susceptibility  maps  were  further  generated 
with  a  threshold  of  0.1.  Susceptibility  values  of  the  veins  were  measured  and  their  relative  errors  were 
calculated. 

c)  Results 

The  Laplacians  of  the  phase  obtained  using  the  discrete  and  continuous  operators  are  shown  in  Fig.  29.  The 
major  differences  in  the  Laplacians  as  well  as  in  the  reconstructed  local  phase  images  are  associated  with  the 
veins.  This  difference  is  also  seen  in  Fig.  30,  which  shows  that  the  phase  image  processed  using  the 
continuous  operator  has  errors  near  the  veins.  After  susceptibility  maps  were  generated  using  both  the 
processed  phase  images  (e.g.,  Figs.  29d  and  29e),  susceptibility  values  of  the  veins  (as  shown  in  Fig.  29d  and 
29e,  excluding  the  superior  sagittal  sinus)  were  measured  as  0.41  ±  0.03  ppm  and  0.36  ±  0.04  ppm,  for  phase 
images  processed  using  discrete  and  continuous  operators,  respectively.  Using  the  discrete  operator,  the 
relative  error  in  the  measured  susceptibility  value  is  9%,  compared  with  the  input  value  of  0.45  ppm.  The 
relative  error  is  20%  for  the  continuous  operator. 
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Figure  29:  (a)  Laplacian  of  the  phase  calculated  using  the  discrete  Laplacian  operator,  (b)  Laplacian 
calculated  using  the  continuous  operator,  (c)  The  resulting  image  of  (b)  subtracted  from  (a),  (d)  Local  phase 
images  obtained  using  the  discrete  Laplacian  operator,  (e)  Phase  image  using  the  continuous  operator,  (f)  The 
resulting  image  of  (e)  subtracted  from  (d). 


d  e 

Figure  30:  (a)  The  simulated  phase  image,  (b)  Phase  image  generated  using  the  discrete  operator  (c)  The 
resulting  image  of  (b)  subtracted  from  (a),  (d)  Phase  image  generated  using  the  continuous  operator,  (e)  The 
resulting  image  of  (d)  subtracted  from  (a). 
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d)  Discussion 


The  Laplacian  of  phase  can  be  calculated  using  the  discrete  Laplacian  kernel  and  phase  differences  between 
pixels,  or  using  the  continuous  Laplacian  operator  based  on  sine  and  cosine  functions  and  Fourier  transforms. 
In  recent  studies,  we  see  a  mixture  of  these  two  methods.  For  example,  Schofield  and  Zhu  [2003]  used  the 
continuous  method  for  phase  unwrapping.  They  claim  that  using  the  Fourier  transform  property  to  calculate  the 
derivatives  and  using  the  sine  and  cosine  functions  help  to  avoid  noise  amplifications.  In  this  study,  random 
noise  was  not  considered.  Hence  the  propagation  of  noise  was  not  fully  evaluated.  Since  for  phase 
unwrapping,  it  is  the  number  of  multiples  of  2n  that  need  to  be  determined,  rather  than  the  local  phase 
information,  it  might  not  have  much  error  associated  with  veins,  as  long  as  the  absolute  value  of  Laplacian  of 
the  local  phase  information  does  not  exceed  n  (if  it  exceeds,  the  number  of  multiples  of  2 n  may  not  be 
accurately  determined).  However,  if  the  Laplacian  is  directly  used  in  phase  image  processing,  the  discrete 
Laplacian  operator  should  still  be  used.  Phase  images  processed  using  the  continuous  operator  may  lead  to 
an  additional  10%  error  in  our  estimated  susceptibility.  Given  the  fact  that  MRI  images  are  discrete  signal,  it  is 
natural  and  more  consistent  to  use  the  discrete  forms. 

It  should  also  be  noted  that,  once  the  Laplacian  kernel  is  directly  used  in  SHARP,  it  is  equivalent  to  a  kernel 
size  of  1  pixel.  This  creates  other  problems  related  to  the  deconvolution  step  in  SHARP,  and  will  lead  to  larger 
errors  in  the  processed  phase  images.  Generally  speaking,  SHARP  works  best  with  filter  sizes  of  6-10  pixels, 
as  shown  in  previous  results. 
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KEY  RESEARCH  ACCOMPLISHMENTS 


The  following  list  of  key  research  accomplishments  is  based  on  the  “Statement  of  Work”  described  in  our 

funded  grant.  In  this  first  year, 

•  We  have  simulated  a  numerical/digital  brain  model  through  forward  simulations  with  known  magnetic 
susceptibility  values  assigned  to  1 1  tissues  inside  the  brain.  This  simulation  has  included  the  partial  volume 
effect,  Gibbs  ringing,  and  tissue  spin  densities  in  magnitude  MR  images. 

•  We  have  built  several  4-straw  or  3-straw  agar  gel  phantoms  with  different  concentrations  of  nanoparticles 
and  gadolinium  mixed  with  distilled  water  and/or  gelatin  in  straws.  Different  diameters  of  straws  were  tried. 

•  We  have  imaged  the  nanoparticle  and  gadolinium  phantoms  at  two  different  orientations  with  an  11 -echo 
susceptibility  weighted  imaging  (SWI)  sequence  and  a  14-echo  spin  echo  sequence  on  our  3T  MRI 
machine.  The  two  different  orientations  include  the  cases  when  straws  are  parallel  and  perpendicular  to  the 
main  field. 

•  We  have  built  the  SHARP  program  originally  developed  by  Schweser  et  al.  [2011], 

•  We  have  polished  the  high-pass  filter  program  and  the  quadratic  fit  program  for  the  background  phase 
removal  from  MR  images. 

•  We  have  built  our  own  algorithm  to  remove  the  background  phase. 

•  We  have  found  that  we  can  relatively  easy  to  quantify  susceptibility  from  larger  straws  with  our  CISSCO 
method  applied  to  MR  images.  We  can  also  quantify  susceptibility  with  CISSCO  from  smaller  straws  with 
larger  uncertainties  plus  the  wall  effect.  Thus  we  have  decided  to  use  straws  with  diameters  of  roughly  6.3 
mm  and  to  prepare  solutions  of  each  material  with  the  following  targeted  magnetic  susceptibility  values:  2 
ppm,  1  ppm,  0.5  ppm,  and  0.25  ppm  in  SI  system. 

•  We  have  applied  SHARP,  high-pass  filter,  quadratic  fit,  and  our  own  algorithm  individually  to  SWI  images 
for  removing  the  background  fields.  After  each  of  those  procedures,  we  have  quantified  the  susceptibility  of 
each  straw  from  at  least  4  echoes  from  SWI.  We  have  compared  results  between  those  methods. 

•  We  have  measured  the  concentrations  of  the  aforementioned  nanoparticle  and  gadolinium  solutions  with 
ICP-AES  (inductively  coupled  plasma  atomic  emission  spectroscopy). 

•  We  have  quantified  relaxation  rates  R2  and  R2*  of  each  aforementioned  concentration  from  MR  images 
and  have  estimated  its  associated  uncertainty.  We  have  found  that  R2  and  R2*  are  about  the  same  for 
gadolinium  samples,  but  not  the  same  for  nanoparticle  samples.  Furthermore,  R2  and  R2*  values  depend 
on  whether  gadolinium  or  nanoparticles  are  mixed  with  distilled  water  or  with  gelatin.  While  we  believe  that 
these  results  depend  on  the  diffusion  constants  in  water  and  gel,  the  ratios  between  the  changes  of  R2  or 
R2*  of  substances  in  water  and  in  gel  are  consistently  different.  We  thus  conclude  that  R2  and  R2*  are  not 
reliable  measurements  for  magnetic  susceptibility  of  substances.  Nonetheless,  quantified  susceptibility 
values  agree  well  with  each  other  within  uncertainties  for  the  same  concentration  of  samples.  Therefore, 
we  have  postponed  the  quantifications  of  R2  or  R2*  for  ferritin  and  calcium  solutions. 

•  We  have  built  two  4-straw  agar  gel  phantoms  with  different  concentrations  of  ferritin  and  calcium  mixed 
with  distilled  water  in  straws. 

•  We  have  imaged  the  ferritin  and  calcium  phantoms  at  two  different  orientations  with  an  1 1 -echo  SWI  and  a 
14-echo  spin  echo  sequence  on  our  3T  MRI  machine.  The  two  different  orientations  include  the  cases 
when  straws  are  parallel  and  perpendicular  to  the  main  field. 

•  We  have  quantified  the  susceptibility  of  each  solution  with  CISSCO  from  at  least  4  echo  times,  based  on 
phase  distributions  outside  each  straw,  with  the  applications  of  the  high-pass  filter  and  SHARP  individually 
to  remove  the  background  field. 

•  We  have  also  quantified  the  phase  value  inside  each  straw  to  examine  whether  these  phase  values  agree 
with  the  susceptibility  values  quantified  from  our  CISSCO  method. 

•  We  have  found  that  phase  values  inside  straws  containing  gadolinium,  calcium,  and  ferritin  are  sometimes 
10%  deviated  from  the  phase  values  calculated  from  the  susceptibility  values  obtained  from  CISSCO. 

•  We  have  identified  some  basic  issues  in  the  SHARP  method  and  have  estimated  its  errors  on  our  digital 
brain  model.  In  fact,  we  have  studied  errors  of  quantified  susceptibility  values  with  SHARP,  high-pass  filter, 
and  a  variable  high-pass  filter  on  our  numerical  brain  model. 

•  We  continue  to  improve  different  methods/programs  mentioned  above,  to  develop  our  own  susceptibility 
quantification  method,  and  to  compare  results  from  different  methods. 
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REPORTABLE  OUTCOMES 


Manuscripts,  abstracts,  or  presentations  with  partial  supports  from  this  grant: 

•  W.  Zheng,  H.  Nichol,  S.  Liu,  Y.-C.  N.  Cheng,  and  E.  M.  Haacke,  Measuring  iron  in  the  brain  using 
quantitative  susceptibility  mapping  and  X-ray  fluorescence  imaging,  Neuroimage,  vol.  78,  pp.  68-74,  2013. 

•  Applications  of  Short  Echo  QSM,  a  talk  presented  by  E.  M.  Haacke  in  the  Second  QSM  Workshop,  Cornell 
University,  Ithaca,  New  York,  USA,  July  26,  2013. 

•  CISSCO  Method  for  Measuring  Susceptibility,  a  talk  presented  by  Y.-C.  N.  Cheng  in  the  Second  QSM 
Workshop,  Cornell  University,  Ithaca,  New  York,  USA,  July  27,  2013. 

Funding  applied  for  based  on  work  supported  by  this  award: 

•  GE/NFL,  Title:  Accurate  Magnetic  Susceptibility  and  Volume  Measurements  with  Error  Analyses  of  Each 
Cerebral  Microbleed  and  Tissue  from  MRI  as  Biomarkers  in  the  Mild  Traumatic  Brain. 

Employment  or  research  opportunities  applied  for  and/or  received  based  on  experience/training  supported  by 

this  award: 

•  A  graduate  student  in  Physics,  He  Xie,  has  passed  his  Ph.D.  prospectus  examine  in  early  September  2013 
with  some  of  the  work  described  above.  He  is  now  a  Ph.D.  candidate. 


CONCLUSION 

We  have  accomplished  most  tasks  provided  in  the  Statement  of  Work.  The  above  list  of  key  accomplishments 
can  already  serve  as  the  conclusion.  In  brief,  our  CISSCO  method  can  accurately  quantify  the  magnetic 
moment  of  a  cylindrical  object.  When  the  cross  section  of  the  object  is  known,  we  can  calculate  the  magnetic 
susceptibility.  The  accuracy  achieved  by  CISSCO  on  the  object  with  a  diameter  of  a  few  pixels  cannot  be 
achieved  by  any  other  susceptibility  quantification  methods.  In  addition,  we  have  studied  the  effects  due  to  four 
background  phase  removal  methods.  Currently  we  have  found  the  inconsistencies  between  phase  values 
inside  a  straw  and  expected  values  estimated  from  quantified  susceptibilities.  We  have  shown  that  the  issues 
can  be  due  to  the  use  of  the  high-pass  filter  and  we  further  suspect  that  the  issues  can  also  be  due  to  the  use 
of  SHARP  (one  of  the  background  phase  removal  method).  We  will  conduct  a  few  more  simulations  and 
resolve  this  issue.  Furthermore,  we  have  found  that  the  relaxation  rate  (R2  or  R2*)  measurements  are  not 
consistent  and  depending  on  the  solid  or  liquid  state  of  a  sample,  the  relaxation  rates  can  be  different,  even  for 
the  same  concentration  of  a  material.  This  will  be  important  for  the  MR  community  to  be  aware.  We  will  also 
need  to  spend  time  to  write  and  publish  our  presented  work  above. 

Our  numerical  brain  model  will  continue  to  be  a  tool  for  tests.  In  the  next  annual  period,  as  listed  in  the 
Statement  of  Work,  we  will  use  or  build  molds  to  prepare  samples  with  arbitrary  geometries  and  quantify  the 
susceptibility  of  each  sample.  We  will  also  need  to  improve  our  own  background  phase  removal  method  to  a 
3D  method.  In  addition  to  existing  susceptibility  quantification  methods,  we  also  want  to  develop  a  better 
method.  These  will  be  challenging  tasks  in  this  coming  annual  period. 
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The  paper  by  Zheng  et  al.  is  attached. 
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Measuring  iron  content  in  the  brain  has  important  implications  for  a  number  of  neurodegenerative  diseases. 
Quantitative  susceptibility  mapping  (QSM),  derived  from  magnetic  resonance  images,  has  been  used  to  mea¬ 
sure  total  iron  content  in  vivo  and  in  post  mortem  brain.  In  this  paper,  we  show  how  magnetic  susceptibility 
from  QSM  correlates  with  total  iron  content  measured  by  X-ray  fluorescence  (XRF)  imaging  and  by  induc¬ 
tively  coupled  plasma  mass  spectrometry  (ICPMS).  The  relationship  between  susceptibility  and  ferritin  iron 
was  estimated  at  1.10  db  0.08  ppb  susceptibility  per  pg  iron/g  wet  tissue,  similar  to  that  of  iron  in  fixed  (fro¬ 
zen/thawed)  cadaveric  brain  and  previously  published  data  from  unfixed  brains.  We  conclude  that  magnetic 
susceptibility  can  provide  a  direct  and  reliable  quantitative  measurement  of  iron  content  and  that  it  can  be 
used  clinically  at  least  in  regions  with  high  iron  content. 

©  2013  Elsevier  Inc.  All  rights  reserved. 


Introduction 

Iron  is  an  important  endogenous  biomarker  for  many  neurological 
diseases  and  normal  aging  (Haacke  et  al.,  2005;  Schenck  and 
Zimmerman,  2004).  Previous  histological  work  has  shown  that  focally 
elevated  iron  deposition  is  associated  with  various  neurological  and 
psychiatric  disorders,  including  multiple  sclerosis  (MS)  (LeVine,  1997), 
Alzheimer's  disease  (Bouras  et  al.,  1997;  Hallgren  and  Sourander, 
1960;  LeVine,  1997),  Huntington's  disease  (Chen  et  al.,  1993;  Dexter  et 
al.,  1991)  and  Parkinson's  disease  (Chen  et  al.,  1993;  Dexter  et  al., 
1991).  Increased  iron  accumulation  has  been  detected  in  chronic  hem¬ 
orrhage,  MS  lesions,  cerebral  infarction,  anemia,  thalassemia,  hemochro¬ 
matosis,  and  NBIA  (neurodegeneration  with  brain  iron  accumulation) 
(Haacke  et  al.,  2005).  An  in  vivo  non-invasive  and  quantitative  estima¬ 
tion  of  non-heme  iron  deposition  (predominantly  ferritin)  is  essential 
to  understand  the  cause  of  iron  accumulation  and  its  distribution  pat¬ 
terns  as  well  as  its  physiological  role  in  any  given  disease  (Bartzokis  et 
al.,  2007;  Gerlach  et  al.,  1994;  Ke  and  Qian,  2003). 

A  variety  of  methods  have  been  used  in  the  past  to  quantify  iron  using 
magnetic  resonance  imaging  (MRI)  (Haacke  et  al.,  2005).  The  standard 
workhorses  in  this  area  are  T2  (House  et  al.,  2007;  Jensen  et  al.,  2010; 
Mitsumori  et  al.,  2012)  and  T2*  (or  R2*  =  1/T2*)  imaging  methods 
that  create  T2*  or  R2*  maps  derived  from  multi-echo  gradient  (recalled) 
echo  magnitude  images.  The  latter  are  particularly  useful  since  gradient 
echo  sequences  are  very  sensitive  to  the  local  susceptibility  induced 
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magnetic  field  inhomogeneity  due  to  iron  (Bartzokis  et  al.,  1993; 
Haacke  et  al.,  1989,  2005;  Ordidge  et  al.,  1994;  Peters  et  al.,  2007; 
Reichenbach  et  al.,  1997).  Further,  T2*  or  R2*  maps  provide  an  important 
contrast  mechanism  to  investigate  brain  tissue  microstructure  and  to 
detect  abnormal  levels  of  brain  iron  (Bartzokis  et  al.,  2007;  Bouras  et  al., 
1997;  Chen  et  al.,  1993;  Dexter  et  al.,  1991;  Haacke  et  al.,  2005,  2009; 
Hallgren  and  Sourander,  1960;  LeVine,  1997;  Wallis  et  al.,  2008). 

In  this  paper,  we  focus  on  susceptibility  measurements  from  phase 
images.  Phase  has  been  used  as  a  means  to  measure  iron  content 
(Haacke  et  al.,  2007).  However,  phase  is  dependent  on  the  geometry 
of  the  object  and  so  it  can  be  misinterpreted.  The  solution  lies  in  using 
a  susceptibility  map  reconstructed  from  the  phase  information.  In 
theory,  this  approach  is  independent  of  field  strength,  echo  time,  the 
object's  relative  orientation  to  the  main  field  and  the  object's  shape 
(Cheng  et  al.,  2009b;  de  Rochefort  et  al.,  2010;  Haacke  et  al.,  2010; 
Kressler  et  al.,  2010;  Li  et  al.,  2011;  Liu  et  al.,  2009;  Marques  and 
Bowtell,  2005;  Schweser  et  al.,  2011;  Shmueli  et  al.,  2009;  Wharton 
and  Bowtell,  2010;  Yao  et  al.,  2009).  Recent  work  has  suggested  that 
susceptibility  changes  in  the  basal  ganglia,  thalamus  and  other  deep 
gray  matter  nuclei  have  better  correlation  with  iron  concentration 
than  phase  information  (Bilgic  et  al.,  2012;  Fukunaga  et  al.,  2010; 
Langkammer  et  al.,  2012b;  Schweser  et  al.,  2011,  2012;  Shmueli  et  al., 
2009;  Wharton  and  Bowtell,  2010;  Yao  et  al.,  2009)  and,  therefore, 
quantitative  susceptibility  mapping  (QSM)  may  provide  a  good  means 
to  study  tissue  iron  content. 

Currently,  the  neuroscience  community  relies  upon  the  50  year  old 
data  on  iron  in  cadaveric  brains  published  by  Hallgren  and  Sourander 
(Hallgren  and  Sourander,  1958).  Total  iron  in  cadaveric  brain  has  been 
measured  using  synchrotron  X-ray  fluorescence  (XRF)  iron  mapping 
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Table  1 

Methodology  and  data  processing. 
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(Haacke  et  al.,  2010) 

(Hopp  et  alM  2010;  Zheng  et  al.,  2012),  proton-induced  X-ray  emission 
mapping  (Butz  et  al.,  2000),  inductively  coupled  plasma  mass  spec¬ 
trometry  (ICPMS)  measurements  (Langkammer  et  al.,  2010,  2012a) 
and  atomic  absorption  spectrometry  measurements  (House  et  al., 
2007).  Among  these,  the  first  two  techniques  can  provide  a  voxel 
by  voxel  quantification  of  iron  content  which  can  then  be  compared 
with  MR  iron  quantification. 

In  this  paper,  our  goal  is  to  develop  an  absolute  quantification 
scale  by  separating  the  iron  induced  susceptibility  change  from 
other  potential  sources  by  comparing  ferritin-gelatin  phantoms  with 
quantified  XRF  iron  maps  of  basal  ganglia  from  cadaver  brains  and 
ICPMS  iron  values. 

Materials  and  methods 

Preparation  of  ferritin  phantoms 

Horse  spleen  ferritin  (Ref.  F4503,  Sigma-Aldrich,  USA)  was  used  to 
prepare  ferritin-gelatin  phantoms.  The  iron  concentration  as  deter¬ 
mined  by  the  supplier  using  ICPMS  was  7.13  ±  0.15  mg/ml.  The  fer¬ 
ritin  solution  was  first  diluted  by  adding  4  ml  of  original  solution 
with  16ml  warm  7%  gelatin  resulting  in  a  stock  solution  with  iron 
concentration  of  about  1426  ±  30  jug/ml.  This  stock  solution  was  se¬ 
rially  diluted  six  times  in  warm  gelatin  by  a  factor  of  2  each  time. 
The  ferritin-doped  gelatin  solutions  as  well  as  pure  gelatin  were  load¬ 
ed  into  straws  and  then  embedded  in  a  pure  gelatin  matrix.  Total  iron 
was  measured  in  aliquots  of  the  ferritin-doped  gelatin  by  XRF  and 
ICPMS.  See  the  detailed  scheme  of  the  experiment  in  Table  1. 

Rapid  scanning  X-ray  fluorescence  (RS-XRF) 

All  XRF  measurements  were  conducted  at  the  Stanford  Synchro¬ 
tron  Radiation  Lightsource  (SSRL).  RS-XRF  images  of  ferritin  phan¬ 
toms  and  cadaveric  brain  were  acquired  at  wiggler  beam  line  10-2 
at  SSRL.  The  samples  were  mounted  onto  a  set  of  motorized  stages 
oriented  at  45°  to  the  incident  beam.  The  incident  beam  (12  keV) 
passing  through  a  tantalum  aperture  produced  a  100  pm  x  100  pm 
spot  on  the  sample  which  was  raster-scanned  in  the  beam  using  a 
dwell  time  of  15  ms/point.  Fluorescent  energy  windows  were  cen¬ 
tered  for  Fe  (6.21-6.70  keV)  as  well  as  all  other  biologically  interest¬ 
ing  elements,  scatter  and  total  incoming  counts.  Elements  were 
quantified  in  pg  iron/g  wet  tissue  by  comparison  of  signal  strength 
with  XRF  calibration  standards  (±5%  uncertainty)  (Micromatter, 
Vancouver,  BC,  Canada)  according  to  Hopp  and  colleagues  (Hopp  et 
al.,  2010)  using  Sam's  Microanalysis  kit  (Webb,  2010).  An  area  of  the 
ferritin-doped  gelatin  block  was  mapped  and  average  counts  were 
compared  with  XRF  calibration  standards. 

Inductively  coupled  plasma  mass  spectrometry 

To  confirm  the  total  iron  content  of  the  ferritin  phantoms,  5  ml 
samples  were  taken  from  the  straws  after  MR  imaging  and  the  iron 
content  was  determined  by  ICPMS  using  an  ELAN  9000  system 
(PerkinElmer,  Waltham,  MA,  USA)  (American  Environmental  Testing 


Fig.  1.  Photograph  of  the  cadaveric  brain  sample  in  gelatin. 


Laboratory  Inc.,  California).  The  samples  were  diluted  to  the  range  ac¬ 
ceptable  for  ICPMS  via  serial  dilutions. 

Preparation  of  the  cadaveric  brain  sample 

One  frozen  coronal  section  (96  mm  long  x  132  mm  wide  x  5  mm 
thick)  of  human  cadaveric  multiple  sclerosis  (MS)  brain  (MS  3852) 
(see  Fig.  1)  was  obtained  from  the  Human  Brain  and  Spinal  Fluid  Re¬ 
source  Center,  Los  Angeles,  CA,  under  the  University  of  Saskatchewan 
ethics  approval  BioREB  06-250.  Coronal  sections  showed  extensive  ir¬ 
regular  demyelination  throughout  the  brainstem.  There  were  also  a 
few  small  scattered  demyelinating  periventricular  foci  (bilateral).  The 
surface  of  the  sample  (a  5  mm  thick  section)  showed  patchy  areas  of 
slight  rarefaction  without  significant  axonal  loss  or  change  in  oligoden¬ 
drocyte  density.  There  were  varying  degrees  of  associated  gliosis.  The 
areas  of  rarefaction  were  associated  with  extensive  demyelination.  To 
reduce  storage  artifacts  such  as  leaching  of  metals,  fresh  autopsy  brain 
was  flash  frozen  and  the  slices  were  shipped  on  dry  ice  and  stored  fro¬ 
zen  until  they  were  thawed  by  immersion  in  buffered  formalin.  After 
6  h  of  fixation,  the  brain  slice  was  drained  and  sealed  in  plastic  prior 
to  initial  synchrotron  imaging  of  the  surface  of  the  slice.  To  resolve  re¬ 
gions  of  interest,  the  slices  were  embedded  in  gelatin  for  MR  imaging. 
The  brain  hemispheres  were  sectioned  to  expose  the  region  of  interest 
and  then  the  slice  was  sealed  in  metal-free  thin  polypropylene  film. 
RS-XRF  images  were  acquired  and  quantified  at  SSRL  (see  the  detailed 
scheme  of  the  experiment  in  Table  1 ). 

MR  imaging  and  image  processing 

Imaging  and  phase  processing  of  ferritin  samples 

MR  data  of  ferritin  samples  were  collected  on  a  3  T  Siemens  Verio 
system  using  a  multi-echo  susceptibility  weighted  imaging  (SWI)  se¬ 
quence  with  11  echoes  (TR  =  40  ms,  FA  =  15°).  The  resolution  was 
1  mm  x  1  mm  x  1  mm  with  a  matrix  of  256  x  256  x  128.  The 
shortest  echo  time  was  5  ms  with  a  2.39  ms  increment  for  the  other 
10  echoes.  Magnitude  and  phase  images  were  reconstructed  from 
the  raw  data  for  each  individual  and  combined  channel.  The  geometry 
of  the  ferritin  samples  was  segmented  from  multi-echo  spin  echo  im¬ 
ages  (TR  =  2000  ms,  resolution  0.22  mm  x  0.22  mm  x  3  mm). 

In  order  to  reconstruct  a  susceptibility  map,  a  pristine  phase  map 
was  required.  That  is,  the  phase  was  unwrapped  and  all  spurious 
phase  information  was  removed.  Phase  images  (TE  =  21.73  ms) 
were  unwrapped  using  Prelude  in  FSL  (Jenkinson,  2003).  To  remove 
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Fig.  2.  Removing  the  background  phase  (TE  =  21.73  ms).  A)  Geometry  of  the  straws  segmented  from  the  spin  echo  images.  B)  Original  phase.  C)  Background  phase  after  extrap¬ 
olation  of  magnetic  fields  into  the  straw  regions.  D)  Subtraction  of  C  from  B  to  reveal  pristine  dipole  effects  due  to  the  iron  in  the  straws. 


the  low  spatial  frequency  background  field  effects,  phase  from  regions 
outside  the  straws  were  chosen,  where  there  were  minimal  remnant  di¬ 
pole  effects.  First,  a  circular  mask  with  a  radius  three  times  that  of  the 
straw  was  defined  and  centered  on  each  straw  and  all  the  information 
inside  this  mask  was  removed  from  the  images.  The  remaining  signal 
was  fit  with  a  quadratic  function  and  extrapolated  back  into  the  masked 
region.  Then  the  estimated  dipole  phase  was  obtained  by  subtracting 
this  modified  background  phase  from  the  original  phase.  The  suscepti¬ 
bility  inside  each  of  the  ferritin  straws  was  assumed  to  be  uniform 
and  was  estimated  using  a  least  squares  fitting  of  the  forward  simulated 
dipole  phase  with  the  estimated  phase  (Neelavalli  et  al.,  2009).  All  the 
steps  were  performed  in  MATLAB  R2009a.  The  results  of  each  step  are 
shown  in  Fig.  2. 

Imaging  and  image  processing  of  cadaveric  brain 

MR  images  were  collected  on  a  3  T  Siemens  Verio  system  using  the 
same  1 1  echo  SWI  sequence  but  with  different  imaging  parameters.  The 
coronal  images  were  acquired  with  a  resolution  0.5  mm  x  0.5  mm  in 
phase  encoding  and  readout  direction  and  0.7  mm  in  the  slice  select 
direction  with  a  readout  bandwidth  of  465  Hz/pixel,  a  field-of-view 
of  256  mm  x  192  mm  with  Nx  =  512,  Ny  =  384  and  Nz  =  40.  The 
shortest  echo  time  was  5.68  ms  with  a  2.57  ms  increment  for  the 
other  10  echoes.  MR  phase  images  (TE  =  8.25  ms)  were  first 
unwrapped  using  Prelude  in  FSL  (Jenkinson,  2003)  and  then  the  back¬ 
ground  phase  was  removed  using  TSVD-SHARP  (Schweser  et  al., 
201 1 )  with  a  kernel  size  of  5  mm.  An  initial  estimation  of  the  suscepti¬ 
bility  distribution  was  obtained  using  truncated  k-space  division,  with  a 
threshold  value  of  0.1.  Due  to  the  presence  of  some  air  bubbles  near  the 
brain  tissue,  the  streaking  artifacts  would  mask  several  important  re¬ 
gions  in  the  susceptibility  map.  Thus,  the  air  bubbles  were  first 
extracted  from  the  susceptiblity  map  by  setting  a  threshold,  since  air 
has  a  much  higher  susceptibility  relative  to  water  than  that  of  brain 
tissue.  The  extracted  susceptibility  maps  of  the  air  bubbles  were  used 
to  predict  their  induced  field  variation  through  a  forward  field  calcula¬ 
tion.  Finally,  the  predicted  fields  induced  by  the  air  bubbles  were  re¬ 
moved  from  the  SHARP  (Schweser  et  al.,  2012)  processed  field  map. 
The  central  region  of  these  air  bubbles  in  phase  images  was  set  to  be 
zero,  in  order  to  reduce  the  streaking  artifacts  caused  by  the  noise  inside 
the  bubble.  This  newly  processed  field  map  was  used  to  generate  the 


final  susceptibility  maps,  using  a  truncated  k-space  division  with  a 
threshold  of  0.1  (Haacke  et  al.,  2010)  via  SPIN  (Signal  Processing  in 
NMR,  Detroit,  MI,  USA)  software. 

Results 

Correlation  between  susceptiblity  and  ferritin  iron  content 

The  susceptibilities  (TE  =  21.73  ms)  of  the  five  empty  straws 
embedded  in  gelatin  were  estimated  at  (9.46  ±  0.015;  9.64  ± 
0.015;  9.46  ±  0.016;  9.65  ±  0.013;  9.46  ±  0.015)  ppm.  Assuming 
that  the  susceptibility  difference  between  the  air  and  gel  is  9.4  ppm 
(Cheng  et  al.,  2009a),  the  total  susceptiblity  measurement  including 
the  background  removal,  straw  geometry  segmentation  error  and 
least  squares  fitting  had  a  bias  of  1.42%. 

The  measured  susceptibilities  (TE  =  21.73  ms)  and  iron  concentra¬ 
tions  of  the  six  ferritin  samples  are  listed  in  Table  2.  The  dipolar  phase 
pattern  outside  the  straw  from  the  sample  with  the  lowest  iron  concen¬ 
tration  (39  ±  6  jig  Fe/ml)  had  its  sign  reversed  compared  with  other 
samples.  This  sample  shows  a  negative  susceptibility  of  —  14ppb 
when  using  the  forward  fitting  approach.  One  possible  explanation  for 
this  could  be  a  small  baseline  shift  coming  from  an  imperfect  back¬ 
ground  removal.  Since  the  iron  concentration  range  that  can  be  mea¬ 
sured  with  XRF  is  broad,  there  was  no  need  for  dilution.  In  contrast, 
ICPMS  requires  dilution  of  samples  to  make  iron  concentration  in  the 
proper  range  for  analysis.  The  results  in  Table  2  show  that  the  iron  con¬ 
tent  measured  by  two  approaches  (XRF  and  ICPMS)  was  essentially  the 
same.  The  correlation  slopes  in  Fig.  3  obtained  from  ICPMS  (1.11  ± 
0.06  ppb  per  jug  iron/ml)  and  XRF  imaging  (1.10  ±  0.08  ppb  per  jug 
iron/ml)  were  close  and  both  were  less  than  the  theoretical  estimation 
of  1.27ppb  per  jug  iron/ml  from  Schenck  (1992). 

Correlation  between  susceptibility  and  iron  in  cadaveric  brain 

In  order  to  correlate  the  susceptibility  and  XRF  iron  maps,  images 
from  both  methods  were  co-registered  (Fig.  4).  ROIs  marked  in  each 
image  were  used  for  a  voxel  by  voxel  comparison  of  susceptibility 
and  iron  measurements  (Table  3).  At  TE  =  8.25  ms,  the  correla¬ 
tion  equations  were  found  to  be  Y  =  0.80(±0.01)  (ppb  susceptibility 


Table  2 

Susceptibilities  of  ferritin  phantoms  as  quantified  from  SWI  phase  data  (TE  =  21.73  ms)  and  iron  concentrations  measured  by  XRF  and  ICPMS.  Data  are  shown  as  mean  ±  one  stan¬ 
dard  deviation. 


Sample  no.  1 

Sample  no.  2 

Sample  no.  3 

Sample  no.  4 

Sample  no.  5 

Sample  no.  6 

Gelatin  solution  (7%) 

Susceptibility  (ppb)  (N  =  19,205) 

XRF  iron  concentration  (pg/ml)  (N  =  961) 
ICPMS  iron  concentration  (pg/ml) 

840  ±  2.4 

790  ±  94 

772  ±  115 

428  ±  1.3 

395  ±  44 

448  ±  67 

271  ±  0.9 

229  ±  32 

240  ±  36 

101  ±  0.4 

110  ±  27 

127  ±  19 

39  ±  0.3 

77  ±  16 

66  ±  10 

-14  ±  0.2 
(Not  available) 
39  ±  6 

N. A.  for  forward  fitting 
(Not  available) 

O. 23  ±  0.11 

Standard  deviation  includes  the  spatial  distribution  variation  in  the  straws. 
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Fig.  3.  Correlation  between  susceptibility  measured  by  MRI  and  total  iron  measured  by 
ICPMS  and  XRF  for  ferritin  phantoms. 


per  jug  iron/g  wet  tissue)  *  X  (jug  iron/g  wet  tissue)  +  10.87(±2.9) 
(ppb  susceptibility)  and  Y  =  0.79(±0.02)  *  X  —  3.66(±4.2)  (ppb 
suscep  rowsep=MlMtibility)  for  left  and  right  hemisphere,  respectively, 
as  shown  in  Fig.  5  (A,  B).  The  phase  images  atTE  =  21.1  ms  were  also 


processed,  the  fitted  equations  were  Y  =  0.78(±0.02)  *  X  —  4.36(±4) 
(ppb  susceptibility)  and  Y  =  0.79(±0.01)  *  X  —  5.22(±2.8)  (ppb  sus¬ 
ceptibility)  for  left  and  right  hemispheres  respectively.  The  slopes  (0.80 
and  0.79  ppb  susceptibility  per  pg  iron/g  wet  tissue)  determined  from 
the  TE  =  21.1  ms  data  were  similar  to  those  from  TE  =  8.25  ms  (0.78 
and  0.79).  Although  phase  is  clearly  modified  as  a  function  of  echo  time, 
tissue  susceptibility  change  is  expected  to  be  and  here  is  shown  to  be 
independent  of  echo  time  (Haacke  et  al.,  2010).  The  estimated  susceptibil¬ 
ity  based  on  our  simulation  of  the  inverse  process  using  the  structures  of  a 
similar  size  showed  an  underestimate  or  bias  of  — 14%. 

Discussion 

Using  ferritin  phantoms  and  a  cadaveric  brain  sample,  we  have 
found  that  the  susceptibility  correlates  reasonably  well  with  the  iron 
measured  by  XRF  and/or  ICPMS  (Fig.  3).  The  cadaveric  brain  used  in 
the  study  was  from  a  person  with  multiple  sclerosis.  It  is  commonly  as¬ 
sumed  that  the  iron  in  normal  and  pathological  MS  brains  is  predomi¬ 
nantly  stored  in  the  form  of  ferritin.  As  long  as  this  assumption  holds, 
the  MS  pathology  will  not  affect  the  susceptibility/iron  correlation 
slope.  Our  correlation  of  iron  content  with  susceptibility  for  cadaveric 
brains  (Fig.  5)  was  comparable  with  that  obtained  by  Langkammer 
et  al.  (2012b).  This  is  expected  since  we  used  SHARP  with  the  same 
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Fig.  4.  Iron  quantified  from  XRF  Fe  mapping  (A,  B)  for  left  and  right  hemispheres;  putative  iron  quantified  as  susceptibility  (TE  =  8.25  ms)  (C,  D).  Images  are  co-registered  and  the 
ROIs  used  for  a  pixel  by  pixel  correlation  are  outlined  in  both  images.  CN:  caudate  nucleus.  PUT:  putamen.  GP:  globus  pallidus.  ROIs  were  defined  by  excluding  the  edges  in  the  map 
for  each  structure. 
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Table  3 

Average  susceptibility  of  a  cadaveric  brain  as  quantified  from  SWI  phase  data  (TE  =  8.25  ms)  and  Fe  measured  using  XRF  imaging.  ROIs  are  defined  in  Fig.  4.  Data  are  shown  as 
mean  ±  one  standard  deviation. 


CN  (left) 

PUT  (left) 

GP  (left) 

CN  (right) 

PUT  (right) 

GP  (right) 

Susceptibility  (ppb) 

Iron  estimated  by  XRF  (pg/g  wet  tissue) 

111  ±  25 

153  ±  28 

152  ±  18 

210  ±  35 

273  ±  73 

338  ±  73 

105  ±  31 

135  ±  28 

121  ±  24 

179  ±  26 

242  ±  40 
277  ±  40 

Mean  is  estimated  as  the  average  within  the  ROI.  Standard  deviation  includes  the  spatial  distribution  variation  within  the  ROI. 


parameters  to  remove  the  background  fields.  The  SWIM  approach  used 
in  this  paper  underestimates  the  susceptibility  by  14%  for  deep  gray 
matter  structures  according  to  our  simulations.  The  homogeneity- 
enabled  incremental  dipole  inversion  (HEIDI)  method  used  by 
Langkammer  et  al.  (2012b)  underestimates  the  susceptibility  by 
about  7%  (Langkammer  et  al.,  2012b;  Schweser  et  al.,  2012).  Our 
slope  (0.8  /  (1-14%)  =  0.93)  is  close  to  that  in  Langkammer  et  al. 
(2012b)  (0.89  /  (1-7%)  =  0.957)  for  deep  gray  matter  when  these  biases 
are  accounted  for.  Since  the  cadaveric  brain  in  our  experiment  was  for¬ 
malin  fixed  and  those  in  Langkammer  et  al.  (2012b)  were  unfixed,  this 
suggests  that  fixation  may  not  change  tissue  susceptibility  in  deep  gray 
matter. 
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Fig.  5.  Correlation  between  susceptibility  and  XRF  iron  measurements  for  all  data 
points  taken  from  each  of  the  regions  demarcated  in  Fig.  4.  A:  fitting  for  left  hemi¬ 
sphere;  B:  fitting  for  right  hemisphere. 


However,  the  slope  of  0.59  ppb  susceptibility  per  pg  iron/g  wet 
tissue  obtained  from  our  in  vivo  data  (Haacke,  2012)  and  other  single 
orientation  results  that  used  Hallgren  and  Sourander's  equation  as  the 
iron  baseline  (Shmueli  et  al.,  2009;  Wharton  and  Bowtell,  2010) 
was  smaller  than  the  0.8  ppb  susceptibility  per  pg  iron/g  wet  tissue 
obtained  from  our  cadaveric  brain  data,  even  though  they  were 
processed  with  the  same  methods.  Thus,  there  appears  to  be  a  difference 
between  in  vivo  and  ex  vivo  susceptibilities  and  their  correlation  with 
iron.  The  reason  for  this  is  unclear  but  could  be  due  to  the  freezing  and 
fixation  process  which  could  affect  local  susceptibilities  of  the  tissue. 

Formalin  fixation  might  change  MR  signal  but  previous  work  on 
myelin  susceptibility  (Lee  et  al.,  2012;  Liu  et  al.,  2011)  demonstrated 
that  the  effect  of  formalin  fixation  on  the  susceptibility  changes  due  to 
myelin  was  subtle.  The  similar  iron/susceptibility  slopes  of  fixed  brain 
in  our  work  and  of  the  unfixed  brains  in  the  work  of  Langkammer  and 
colleagues  (Langkammer  et  al.,  2012b)  further  supports  the  view  that 
formalin  fixation  has  negligible  effect  on  susceptibility.  The  effects  of 
fixation  on  R2  and  thus  R2*  values,  however,  are  known  to  be  substan¬ 
tial  (Dawe  et  al.,  2009;  Lee  et  al.,  2012;  Pfefferbaum  et  al.,  2004; 
Schmierer  et  al.,  2008)  and  are  beyond  the  scope  of  this  paper. 

The  susceptibility/iron  correlation  slopes  obtained  from  cadaveric 
and  in  vivo  brains  in  Table  4  are  generally  smaller  than  the  theoretical 
slope  of  1.32  ppb  susceptibility  per  pg  iron/g  wet  tissue  except  for  the 
slope  obtained  with  myelin  correction  in  Schweser  et  al.  (2011).  One 
possible  reason  for  the  smaller  slope  from  the  in  vivo  human  brains  is 
that  there  are  still  some  forms  of  iron  that  are  MR  invisible  although 
these  may  be  in  other  species  that  are  known  to  be  present  at  low  levels 
(Hopp  et  al.,  2010).  A  second  explanation  for  smaller  slopes  seen  in  our 
work  (Haacke,  2012)  and  other  work  (Shmueli  et  al.,  2009;  Wharton 
and  Bowtell,  2010)  is  that  Hallgren  and  Sourander's  measurements  of 
total  iron  (Hallgren  and  Sourander,  1958)  may  not  be  accurate.  Third, 
susceptibility  mapping  is  known  to  have  a  bias  and  leads  to  a  smaller 
slope,  but  this  bias  can  be  potentially  corrected  (J.  Liu  et  al.,  2012;  T. 
Liu  et  al.,  2012;  Schweser  et  al.,  2012;  Wharton  and  Bowtell,  2010). 
Other  possible  factors  that  have  been  explored  include  contributions 
of  myelin  (Duyn  et  al.,  2007;  Liu  et  al.,  201 1 ;  Ogg  et  al.,  1999),  chemical 
exchange  between  water  and  macromolecular  protons  (Luo  et  al.,  2010; 
Shmueli  et  al.,  201 1 ;  Zhong  et  al.,  2008)  and  microstructure  orientation 
(He  and  Yablonskiy,  2009;  Lee  et  al.,  2010;  Liu,  2010).  Indeed,  it  could 
well  be  a  combination  of  all  these  sources  that  lead  to  different 
measurements  of  iron  in  vivo  and  ex  vivo.  Despite  these  imperfections, 
the  slopes  for  susceptibility  versus  iron  content  are  generally  consistent 
between  both  ex  vivo  studies  (this  paper  and  Langkammer  et  al.,  2012b) 
and  in  vivo  studies  using  similar  susceptibility  mapping  methods 
(see  Table  4)  (Haacke,  2012;  Shmueli  et  al.,  2009;  Wharton  and  Bowtell, 
2010). 

Conclusion 

Our  results  suggest  that  susceptibility  changes  from  iron  measured 
in  ex  vivo  studies  reasonably  reflect  iron  content  even  for  in  vivo  stud¬ 
ies,  although  the  predicted  values  may  be  underestimated.  Our  study 
further  demonstrates  that  the  correlation  of  susceptibility  with  iron  is 
consistent  with  other  results  in  the  literature  and  is  independent  of 
echo  time  and  orientation.  Thus,  susceptibility  would  appear  to  be  a 
direct  and  reliable  quantitative  indication  of  iron,  especially  for  brain 
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Table  4 

Correlation  between  susceptibility  mapping  and  iron  concentration. 


Authors 

Correlation 

slope3 

Structures 

Background  removal 

QSM  method13 

Myelin 

correction 

Field 

(Tesla) 

Sample 

Iron 

This  paper 

1.11 

N.A. 

Quadratic  fitting 

Forward  fitting 

N.A. 

3  T 

Ferritin 

ICPMS 

This  paper 

1.10 

N.A. 

Quadratic  fitting 

Forward  fitting 

N.A. 

3  T 

Ferritin 

XRF 

This  paper 

0.80 

GP,  PUT,  CN 

SHARP 

TKD1 (SO) 

No 

3  T 

MS  cadaveric 
brain  (fixed) 

XRF 

Haacke  (2012) 

0.59 

GP,  PUT,  CN 

SHARP 

TKD1 (SO) 

No 

3  T 

In  vivo  brains 

H&SC 

Shmueli  et  al.  (2009) 

0.56 

PUT,  RN,  SN 

Polynomial  fitting 

TKD2 (SO) 

No 

7T 

In  vivo  brain 

H&S 

Wharton  and  Bowtell  (2010) 

0.75/0.6 

GP,  SN,  RN,  PUT, 
CN,  TH,  GM 

Simulated  geometric 
effect  +  fitting 

TKD3  (MO/SO) 

No 

7T 

In  vivo  brains 

H&S 

Langkammer  et  al.  (2012a) 

0.89 

GP,  PUT,  CN,  TH 

SHARP 

HEIDI  (SO) 

No 

3  T 

Unfixed  cadaveric 
brains 

ICPMS 

Schweser  et  al.  (2011) 

1.30 

GP,  SN,  DN,  PUT, 
CN,  TH,  WM,  GM 

SHARP 

MO  regularization 

Yes 

3  T 

In  vivo  brains 

H&S 

a  The  unit  of  the  slope  for  human  brain  is  susceptibility/pg  iron/g  wet  tissue  (p  =  1 .04g/ml  at  36.5  °C) ;  the  unit  of  the  slope  for  the  ferritin  solution  is  ppb  susceptibility/pg  iron/ml 
and  the  corresponding  theoretical  value  is  1.27. 

b  SO:  single  orientation;  MO:  multiple  orientation;  TKD:  thresholded  k-space  division.  TKD1:  Haacke  et  al.,  2010;  TKD2:  Shmueli  et  al.,  2009;  TKD3:  (Wharton  and  Bowtell,  2010). 
c  H&S:  (Hallgren  and  Sourander,  1958). 


regions  with  high  iron  content.  Susceptibility  mapping  provides  a  reli¬ 
able  tool  for  clinical  investigations  of  iron  that  could  be  used  to  study 
changes  in  iron  over  time  or  within  a  given  age-matched  population. 
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